function r = hydro(Ti, TIME_STEPS)
  % First order hydro
  %  eliminating the u-dot-dot terms 
  %  fsolve search method
  close all;
  fig = figure; 
  load('cm', 'mycmap');
  set(fig, 'Colormap', mycmap);

  HYDROLIZE = false;

  GRID_SIZE = 32;
  DX = 0.01;
  DT = 0.001;
  VECTOR_SUBSAMPLING = 16;
  [gridx, gridy] = meshgrid([1:VECTOR_SUBSAMPLING:GRID_SIZE], [1:VECTOR_SUBSAMPLING:GRID_SIZE]);


  VIDEO = false;
  FPS = 30;
  SECONDS_PER_LIGHTYEAR = 15;
  if (VIDEO)
    vt = VideoWriter('v_t.avi', 'Archival');
    open(vt);
  end

  ZERO_HYDRO = false;



  A = sqrt(3)*pi/18 - log(3)/2;

  ux = zeros(GRID_SIZE, GRID_SIZE);
  uy = zeros(GRID_SIZE, GRID_SIZE);
  T  = zeros(GRID_SIZE, GRID_SIZE);
  
  dtT  = zeros(GRID_SIZE, GRID_SIZE); 
  dtux = zeros(GRID_SIZE, GRID_SIZE); 
  dtuy = zeros(GRID_SIZE, GRID_SIZE); 

  dxT = zeros(GRID_SIZE, GRID_SIZE); 
  dyT = zeros(GRID_SIZE, GRID_SIZE); 
  dxux = zeros(GRID_SIZE, GRID_SIZE); 
  dyux = zeros(GRID_SIZE, GRID_SIZE); 
  dxuy = zeros(GRID_SIZE, GRID_SIZE); 
  dyuy = zeros(GRID_SIZE, GRID_SIZE); 

  dxdxT = zeros(GRID_SIZE, GRID_SIZE); 
  dxdyT = zeros(GRID_SIZE, GRID_SIZE); 
  dydyT = zeros(GRID_SIZE, GRID_SIZE); 

  dxdxux = zeros(GRID_SIZE, GRID_SIZE); 
  dxdyux = zeros(GRID_SIZE, GRID_SIZE); 
  dydyux = zeros(GRID_SIZE, GRID_SIZE); 

  dxdxuy = zeros(GRID_SIZE, GRID_SIZE); 
  dxdyuy = zeros(GRID_SIZE, GRID_SIZE); 
  dydyuy = zeros(GRID_SIZE, GRID_SIZE); 



  sumT = 0;
  % initial conditions for vx vy come here
  % vx
  % vy
  % T
  TEMP_OFFSET =0.7;
  cx = 0.5 * GRID_SIZE + 0.5;
  cy = 0.5 * GRID_SIZE + 0.5;
  for x = 1:GRID_SIZE
    for y = 1:GRID_SIZE
       T(x, y) = 0.004 * exp(-((x-cx)^2 + (y-cy)^2)/12);
      % T(x,y) = 0.2*sin(2*pi*x/GRID_SIZE);
    end
  end
  % T(:,:) = Ti(:,:);
  T = T + TEMP_OFFSET;
  for x = 1:GRID_SIZE
     for y = 1:GRID_SIZE
       %T(x,y) = T(x,y) +  0.01 * rand();
       %vx(x,y) = 0.001 * rand();
       %vy(x,y) = 0.001 * rand();
     end
  end

  %vx = 0.000000 * ones(GRID_SIZE, GRID_SIZE);;
  %vy = 0.2 * ones(GRID_SIZE, GRID_SIZE);;
  
  % ****************************************************************
  %  hydrolize the initial conditions
  % ****************************************************************
  if (HYDROLIZE)
    level = 4;
    k = zeros(GRID_SIZE+1, GRID_SIZE+1);  
    [cx, cy] = meshgrid(0:GRID_SIZE-1,0:GRID_SIZE-1);
    [kx, ky] = meshgrid(-GRID_SIZE/2:GRID_SIZE/2);
    for kx = -GRID_SIZE/level:GRID_SIZE/level
      for ky = -GRID_SIZE/level:GRID_SIZE/level
        ft = T .* exp(-2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);
        k(kx+GRID_SIZE/2+1, ky+GRID_SIZE/2+1) = sum(sum(ft));
      end
    end
    T = zeros(GRID_SIZE:GRID_SIZE);
    for kx = -GRID_SIZE/level:GRID_SIZE/level
      for ky = -GRID_SIZE/level:GRID_SIZE/level
        T = T + k(kx+GRID_SIZE/2+1,ky+GRID_SIZE/2+1) * exp(2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);
        % vx = vx + 20*(rand()-0.5)*exp(2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);
        % vy = vy + 20*(rand()-0.5)*exp(2 * pi * i * ( kx * cx + ky * cy ) / GRID_SIZE);

      end
    end
    T = real(T / GRID_SIZE / GRID_SIZE);
    %ux = real(ux / GRID_SIZE / GRID_SIZE);
    %uy = real(uy / GRID_SIZE / GRID_SIZE);
  end



  % ****************************************************************
  %  compute derivative kernels
  % ****************************************************************
  % can upgrade on this -- trefethen pg 6

  d1 = zeros(GRID_SIZE, GRID_SIZE);
  d2 = zeros(GRID_SIZE, GRID_SIZE);
  for x = 1:GRID_SIZE
    for y=1:GRID_SIZE
      s1 = 0;
      s2 = 0;
      for k = -GRID_SIZE/2 : GRID_SIZE/2
        s1 = s1 + 2. * pi * i / GRID_SIZE / GRID_SIZE * k * exp(2 * pi * i * k * (x - y) / GRID_SIZE);
        s2 = s2 + (2. * pi * i / GRID_SIZE * k)^2 / GRID_SIZE * exp(2 * pi * i * k * (x - y) / GRID_SIZE);
      end
      d1(x,y) = real(s1)/DX;
      d2(x,y) = real(s2)/DX/DX;
    end
  end
  h = 2*pi/GRID_SIZE;
  col1 = [0 .5*(-1).^(1:GRID_SIZE-1).*cot((1:GRID_SIZE-1)*h/2)]  * h/DX;
  max(max(d1 - toeplitz(col1, col1([1 GRID_SIZE:-1:2]))))
  col2 = [(-pi^2/(3*h*h)-1/6) .5*(-1).^(0:GRID_SIZE-2).*(csc((1:GRID_SIZE-1)*h/2).^2)]  * h*h/(DX*DX);
  % d2 = toeplitz(col2);
  

  %dxT = dxdx(T);
  %dxT2 = dxdx2(T);
  % image(imagesc(dxT));
  %figure;
  %image(imagesc(dxT2));
  %image((dxT) * 128  + 128); 
  % image((dxT) * 128 + 128); 
  % truesize(fig);
  %fig=figure;
  %truesize(fig);
  %keyboard();
 % image((T-TEMP_OFFSET)*256+128);

  %return

  dxux = dx(ux);
  dyux = dy(ux);
  dxuy = dx(uy);
  dyuy = dy(uy);
  dxT = dx(T);
  dyT = dy(T);

  dxdxux = dxdx(ux);
  dxdyux = dx(dy(ux));
  dydyux = dydy(ux);
  dxdxuy = dxdx(uy);
  dxdyuy = dx(dy(uy));
  dydyuy = dydy(uy);
  dxdxT = dxdx(T);
  dxdyT = dx(dy(T));
  dydyT = dydy(T);


  tic;


  for t = 1:TIME_STEPS
    disp('time step ');
    disp(t);

    % solving non-linear equation of type:
    % L(u, T, u-dot, T-dot) + R(u, T) == 0  -- 3 components
    % each of L, R has 3 components 

    % precompute the R functions
    % R1
    a = -1;
    R1 = (32.*pi.^3.*T.^3.*(ux.*uy.*dyux + dyuy + ux.^2.*dyuy + 2.*uy.^2.*dyuy + dxux + 2.*ux.^2.*dxux + uy.^2.*dxux + ux.*uy.*dxuy))./(9.*sqrt(1 + ux.^2 + uy.^2)) + (1./(9.*(1 + ux.^2 + uy.^2).^(3./2).*(2 + ux.^2 + uy.^2).^2)).*(8.*pi.^2.*(48.*pi.*T.^2.*uy.*dyT + 144.*pi.*T.^2.*ux.^2.*uy.*dyT + 156.*pi.*T.^2.*ux.^4.*uy.*dyT + 72.*pi.*T.^2.*ux.^6.*uy.*dyT + 12.*pi.*T.^2.*ux.^8.*uy.*dyT + 144.*pi.*T.^2.*uy.^3.*dyT + 312.*pi.*T.^2.*ux.^2.*uy.^3.*dyT + 216.*pi.*T.^2.*ux.^4.*uy.^3.*dyT + 48.*pi.*T.^2.*ux.^6.*uy.^3.*dyT + 156.*pi.*T.^2.*uy.^5.*dyT + 216.*pi.*T.^2.*ux.^2.*uy.^5.*dyT + 72.*pi.*T.^2.*ux.^4.*uy.^5.*dyT + 72.*pi.*T.^2.*uy.^7.*dyT + 48.*pi.*T.^2.*ux.^2.*uy.^7.*dyT + 12.*pi.*T.^2.*uy.^9.*dyT + 4.*ux.^2.*dyT.^2 + 8.*ux.^4.*dyT.^2 + 5.*ux.^6.*dyT.^2 + ux.^8.*dyT.^2 + 4.*uy.^2.*dyT.^2 + 18.*ux.^2.*uy.^2.*dyT.^2 + 18.*ux.^4.*uy.^2.*dyT.^2 + 5.*ux.^6.*uy.^2.*dyT.^2 + 10.*uy.^4.*dyT.^2 + 21.*ux.^2.*uy.^4.*dyT.^2 + 9.*ux.^4.*uy.^4.*dyT.^2 + 8.*uy.^6.*dyT.^2 + 7.*ux.^2.*uy.^6.*dyT.^2 + 2.*uy.^8.*dyT.^2 - 4.*T.*ux.^3.*dyT.*dyux - 4.*T.*ux.^5.*dyT.*dyux - T.*ux.^7.*dyT.*dyux - 8.*T.*ux.*uy.^2.*dyT.*dyux - 18.*T.*ux.^3.*uy.^2.*dyT.*dyux - 8.*T.*ux.^5.*uy.^2.*dyT.*dyux - 14.*T.*ux.*uy.^4.*dyT.*dyux - 13.*T.*ux.^3.*uy.^4.*dyT.*dyux - 6.*T.*ux.*uy.^6.*dyT.*dyux - 4.*T.^2.*dyux.^2 - 4.*T.^2.*ux.^2.*dyux.^2 - T.^2.*ux.^4.*dyux.^2 - 12.*T.^2.*uy.^2.*dyux.^2 - 4.*T.^2.*ux.^2.*uy.^2.*dyux.^2 + 3.*T.^2.*ux.^4.*uy.^2.*dyux.^2 - 11.*T.^2.*uy.^4.*dyux.^2 - 3.*T.^2.*uy.^6.*dyux.^2 + 4.*T.*ux.^2.*uy.*dyT.*dyuy + 8.*T.*ux.^4.*uy.*dyT.*dyuy + 3.*T.*ux.^6.*uy.*dyT.*dyuy + 6.*T.*ux.^2.*uy.^3.*dyT.*dyuy + 4.*T.*ux.^4.*uy.^3.*dyT.*dyuy - 2.*T.*uy.^5.*dyT.*dyuy - T.*ux.^2.*uy.^5.*dyT.*dyuy - 2.*T.*uy.^7.*dyT.*dyuy + 4.*T.^2.*ux.*uy.*dyux.*dyuy - 2.*T.^2.*ux.^3.*uy.*dyux.*dyuy - 3.*T.^2.*ux.^5.*uy.*dyux.*dyuy + 14.*T.^2.*ux.*uy.^3.*dyux.*dyuy + 6.*T.^2.*ux.^3.*uy.^3.*dyux.*dyuy + 9.*T.^2.*ux.*uy.^5.*dyux.*dyuy - 4.*T.^2.*dyuy.^2 - 6.*T.^2.*ux.^2.*dyuy.^2 - 2.*T.^2.*ux.^4.*dyuy.^2 - 10.*T.^2.*uy.^2.*dyuy.^2 - 16.*T.^2.*ux.^2.*uy.^2.*dyuy.^2 - 6.*T.^2.*ux.^4.*uy.^2.*dyuy.^2 - 6.*T.^2.*uy.^4.*dyuy.^2 - 6.*T.^2.*ux.^2.*uy.^4.*dyuy.^2 + 4.*T.*ux.^2.*dydyT + 8.*T.*ux.^4.*dydyT + 5.*T.*ux.^6.*dydyT + T.*ux.^8.*dydyT + 4.*T.*uy.^2.*dydyT + 18.*T.*ux.^2.*uy.^2.*dydyT + 18.*T.*ux.^4.*uy.^2.*dydyT + 5.*T.*ux.^6.*uy.^2.*dydyT + 10.*T.*uy.^4.*dydyT + 21.*T.*ux.^2.*uy.^4.*dydyT + 9.*T.*ux.^4.*uy.^4.*dydyT + 8.*T.*uy.^6.*dydyT + 7.*T.*ux.^2.*uy.^6.*dydyT + 2.*T.*uy.^8.*dydyT - 4.*T.^2.*ux.*dydyux - 8.*T.^2.*ux.^3.*dydyux - 5.*T.^2.*ux.^5.*dydyux - T.^2.*ux.^7.*dydyux - 12.*T.^2.*ux.*uy.^2.*dydyux - 16.*T.^2.*ux.^3.*uy.^2.*dydyux - 5.*T.^2.*ux.^5.*uy.^2.*dydyux - 11.*T.^2.*ux.*uy.^4.*dydyux - 7.*T.^2.*ux.^3.*uy.^4.*dydyux - 3.*T.^2.*ux.*uy.^6.*dydyux - 4.*T.^2.*uy.*dydyuy - 6.*T.^2.*ux.^2.*uy.*dydyuy - 2.*T.^2.*ux.^4.*uy.*dydyuy - 10.*T.^2.*uy.^3.*dydyuy - 10.*T.^2.*ux.^2.*uy.^3.*dydyuy - 2.*T.^2.*ux.^4.*uy.^3.*dydyuy - 8.*T.^2.*uy.^5.*dydyuy - 4.*T.^2.*ux.^2.*uy.^5.*dydyuy - 2.*T.^2.*uy.^7.*dydyuy + 48.*pi.*T.^2.*ux.*dxT + 144.*pi.*T.^2.*ux.^3.*dxT + 156.*pi.*T.^2.*ux.^5.*dxT + 72.*pi.*T.^2.*ux.^7.*dxT + 12.*pi.*T.^2.*ux.^9.*dxT + 144.*pi.*T.^2.*ux.*uy.^2.*dxT + 312.*pi.*T.^2.*ux.^3.*uy.^2.*dxT + 216.*pi.*T.^2.*ux.^5.*uy.^2.*dxT + 48.*pi.*T.^2.*ux.^7.*uy.^2.*dxT + 156.*pi.*T.^2.*ux.*uy.^4.*dxT + 216.*pi.*T.^2.*ux.^3.*uy.^4.*dxT + 72.*pi.*T.^2.*ux.^5.*uy.^4.*dxT + 72.*pi.*T.^2.*ux.*uy.^6.*dxT + 48.*pi.*T.^2.*ux.^3.*uy.^6.*dxT + 12.*pi.*T.^2.*ux.*uy.^8.*dxT + 4.*ux.^3.*uy.*dyT.*dxT + 6.*ux.^5.*uy.*dyT.*dxT + 2.*ux.^7.*uy.*dyT.*dxT + 4.*ux.*uy.^3.*dyT.*dxT + 12.*ux.^3.*uy.^3.*dyT.*dxT + 6.*ux.^5.*uy.^3.*dyT.*dxT + 6.*ux.*uy.^5.*dyT.*dxT + 6.*ux.^3.*uy.^5.*dyT.*dxT + 2.*ux.*uy.^7.*dyT.*dxT - 8.*T.*uy.*dyux.*dxT - 18.*T.*ux.^2.*uy.*dyux.*dxT - 17.*T.*ux.^4.*uy.*dyux.*dxT - 6.*T.*ux.^6.*uy.*dyux.*dxT - 14.*T.*uy.^3.*dyux.*dxT - 24.*T.*ux.^2.*uy.^3.*dyux.*dxT - 13.*T.*ux.^4.*uy.^3.*dyux.*dxT - 7.*T.*uy.^5.*dyux.*dxT - 8.*T.*ux.^2.*uy.^5.*dyux.*dxT - T.*uy.^7.*dyux.*dxT + 8.*T.*ux.*dyuy.*dxT + 22.*T.*ux.^3.*dyuy.*dxT + 19.*T.*ux.^5.*dyuy.*dxT + 5.*T.*ux.^7.*dyuy.*dxT + 18.*T.*ux.*uy.^2.*dyuy.*dxT + 28.*T.*ux.^3.*uy.^2.*dyuy.*dxT + 10.*T.*ux.^5.*uy.^2.*dyuy.*dxT + 9.*T.*ux.*uy.^4.*dyuy.*dxT + 5.*T.*ux.^3.*uy.^4.*dyuy.*dxT + 4.*ux.^2.*dxT.^2 + 10.*ux.^4.*dxT.^2 + 8.*ux.^6.*dxT.^2 + 2.*ux.^8.*dxT.^2 + 4.*uy.^2.*dxT.^2 + 18.*ux.^2.*uy.^2.*dxT.^2 + 21.*ux.^4.*uy.^2.*dxT.^2 + 7.*ux.^6.*uy.^2.*dxT.^2 + 8.*uy.^4.*dxT.^2 + 18.*ux.^2.*uy.^4.*dxT.^2 + 9.*ux.^4.*uy.^4.*dxT.^2 + 5.*uy.^6.*dxT.^2 + 5.*ux.^2.*uy.^6.*dxT.^2 + uy.^8.*dxT.^2 + 8.*T.*uy.*dyT.*dxux + 18.*T.*ux.^2.*uy.*dyT.*dxux + 9.*T.*ux.^4.*uy.*dyT.*dxux + 22.*T.*uy.^3.*dyT.*dxux + 28.*T.*ux.^2.*uy.^3.*dyT.*dxux + 5.*T.*ux.^4.*uy.^3.*dyT.*dxux + 19.*T.*uy.^5.*dyT.*dxux + 10.*T.*ux.^2.*uy.^5.*dyT.*dxux + 5.*T.*uy.^7.*dyT.*dxux - 12.*T.^2.*ux.*uy.*dyux.*dxux - 6.*T.^2.*ux.^3.*uy.*dyux.*dxux + 3.*T.^2.*ux.^5.*uy.*dyux.*dxux - 22.*T.^2.*ux.*uy.^3.*dyux.*dxux - 6.*T.^2.*ux.^3.*uy.^3.*dyux.*dxux - 9.*T.^2.*ux.*uy.^5.*dyux.*dxux + 8.*T.^2.*dyuy.*dxux + 16.*T.^2.*ux.^2.*dyuy.*dxux + 8.*T.^2.*ux.^4.*dyuy.*dxux + 16.*T.^2.*uy.^2.*dyuy.*dxux + 32.*T.^2.*ux.^2.*uy.^2.*dyuy.*dxux + 12.*T.^2.*ux.^4.*uy.^2.*dyuy.*dxux + 8.*T.^2.*uy.^4.*dyuy.*dxux + 12.*T.^2.*ux.^2.*uy.^4.*dyuy.*dxux - 2.*T.*ux.^5.*dxT.*dxux - 2.*T.*ux.^7.*dxT.*dxux + 4.*T.*ux.*uy.^2.*dxT.*dxux + 6.*T.*ux.^3.*uy.^2.*dxT.*dxux - T.*ux.^5.*uy.^2.*dxT.*dxux + 8.*T.*ux.*uy.^4.*dxT.*dxux + 4.*T.*ux.^3.*uy.^4.*dxT.*dxux + 3.*T.*ux.*uy.^6.*dxT.*dxux - 4.*T.^2.*dxux.^2 - 10.*T.^2.*ux.^2.*dxux.^2 - 6.*T.^2.*ux.^4.*dxux.^2 - 6.*T.^2.*uy.^2.*dxux.^2 - 16.*T.^2.*ux.^2.*uy.^2.*dxux.^2 - 6.*T.^2.*ux.^4.*uy.^2.*dxux.^2 - 2.*T.^2.*uy.^4.*dxux.^2 - 6.*T.^2.*ux.^2.*uy.^4.*dxux.^2 - 8.*T.*ux.*dyT.*dxuy - 14.*T.*ux.^3.*dyT.*dxuy - 7.*T.*ux.^5.*dyT.*dxuy - T.*ux.^7.*dyT.*dxuy - 18.*T.*ux.*uy.^2.*dyT.*dxuy - 24.*T.*ux.^3.*uy.^2.*dyT.*dxuy - 8.*T.*ux.^5.*uy.^2.*dyT.*dxuy - 17.*T.*ux.*uy.^4.*dyT.*dxuy - 13.*T.*ux.^3.*uy.^4.*dyT.*dxuy - 6.*T.*ux.*uy.^6.*dyT.*dxuy - 8.*T.^2.*dyux.*dxuy - 16.*T.^2.*ux.^2.*dyux.*dxuy - 12.*T.^2.*ux.^4.*dyux.*dxuy - 3.*T.^2.*ux.^6.*dyux.*dxuy - 16.*T.^2.*uy.^2.*dyux.*dxuy - 8.*T.^2.*ux.^2.*uy.^2.*dyux.*dxuy + 3.*T.^2.*ux.^4.*uy.^2.*dyux.*dxuy - 12.*T.^2.*uy.^4.*dyux.*dxuy + 3.*T.^2.*ux.^2.*uy.^4.*dyux.*dxuy - 3.*T.^2.*uy.^6.*dyux.*dxuy - 12.*T.^2.*ux.*uy.*dyuy.*dxuy - 22.*T.^2.*ux.^3.*uy.*dyuy.*dxuy - 9.*T.^2.*ux.^5.*uy.*dyuy.*dxuy - 6.*T.^2.*ux.*uy.^3.*dyuy.*dxuy - 6.*T.^2.*ux.^3.*uy.^3.*dyuy.*dxuy + 3.*T.^2.*ux.*uy.^5.*dyuy.*dxuy - 8.*T.*ux.^2.*uy.*dxT.*dxuy - 14.*T.*ux.^4.*uy.*dxT.*dxuy - 6.*T.*ux.^6.*uy.*dxT.*dxuy - 4.*T.*uy.^3.*dxT.*dxuy - 18.*T.*ux.^2.*uy.^3.*dxT.*dxuy - 13.*T.*ux.^4.*uy.^3.*dxT.*dxuy - 4.*T.*uy.^5.*dxT.*dxuy - 8.*T.*ux.^2.*uy.^5.*dxT.*dxuy - T.*uy.^7.*dxT.*dxuy + 4.*T.^2.*ux.*uy.*dxux.*dxuy + 14.*T.^2.*ux.^3.*uy.*dxux.*dxuy + 9.*T.^2.*ux.^5.*uy.*dxux.*dxuy - 2.*T.^2.*ux.*uy.^3.*dxux.*dxuy + 6.*T.^2.*ux.^3.*uy.^3.*dxux.*dxuy - 3.*T.^2.*ux.*uy.^5.*dxux.*dxuy - 4.*T.^2.*dxuy.^2 - 12.*T.^2.*ux.^2.*dxuy.^2 - 11.*T.^2.*ux.^4.*dxuy.^2 - 3.*T.^2.*ux.^6.*dxuy.^2 - 4.*T.^2.*uy.^2.*dxuy.^2 - 4.*T.^2.*ux.^2.*uy.^2.*dxuy.^2 - T.^2.*uy.^4.*dxuy.^2 + 3.*T.^2.*ux.^2.*uy.^4.*dxuy.^2 + 4.*T.*ux.^3.*uy.*dxdyT + 6.*T.*ux.^5.*uy.*dxdyT + 2.*T.*ux.^7.*uy.*dxdyT + 4.*T.*ux.*uy.^3.*dxdyT + 12.*T.*ux.^3.*uy.^3.*dxdyT + 6.*T.*ux.^5.*uy.^3.*dxdyT + 6.*T.*ux.*uy.^5.*dxdyT + 6.*T.*ux.^3.*uy.^5.*dxdyT + 2.*T.*ux.*uy.^7.*dxdyT - 6.*T.^2.*ux.^2.*uy.*dxdyux - 9.*T.^2.*ux.^4.*uy.*dxdyux - 3.*T.^2.*ux.^6.*uy.*dxdyux + 2.*T.^2.*uy.^3.*dxdyux - 6.*T.^2.*ux.^2.*uy.^3.*dxdyux - 5.*T.^2.*ux.^4.*uy.^3.*dxdyux + 3.*T.^2.*uy.^5.*dxdyux - T.^2.*ux.^2.*uy.^5.*dxdyux + T.^2.*uy.^7.*dxdyux + 2.*T.^2.*ux.^3.*dxdyuy + 3.*T.^2.*ux.^5.*dxdyuy + T.^2.*ux.^7.*dxdyuy - 6.*T.^2.*ux.*uy.^2.*dxdyuy - 6.*T.^2.*ux.^3.*uy.^2.*dxdyuy - T.^2.*ux.^5.*uy.^2.*dxdyuy - 9.*T.^2.*ux.*uy.^4.*dxdyuy - 5.*T.^2.*ux.^3.*uy.^4.*dxdyuy - 3.*T.^2.*ux.*uy.^6.*dxdyuy + 4.*T.*ux.^2.*dxdxT + 10.*T.*ux.^4.*dxdxT + 8.*T.*ux.^6.*dxdxT + 2.*T.*ux.^8.*dxdxT + 4.*T.*uy.^2.*dxdxT + 18.*T.*ux.^2.*uy.^2.*dxdxT + 21.*T.*ux.^4.*uy.^2.*dxdxT + 7.*T.*ux.^6.*uy.^2.*dxdxT + 8.*T.*uy.^4.*dxdxT + 18.*T.*ux.^2.*uy.^4.*dxdxT + 9.*T.*ux.^4.*uy.^4.*dxdxT + 5.*T.*uy.^6.*dxdxT + 5.*T.*ux.^2.*uy.^6.*dxdxT + T.*uy.^8.*dxdxT - 4.*T.^2.*ux.*dxdxux - 10.*T.^2.*ux.^3.*dxdxux - 8.*T.^2.*ux.^5.*dxdxux - 2.*T.^2.*ux.^7.*dxdxux - 6.*T.^2.*ux.*uy.^2.*dxdxux - 10.*T.^2.*ux.^3.*uy.^2.*dxdxux - 4.*T.^2.*ux.^5.*uy.^2.*dxdxux - 2.*T.^2.*ux.*uy.^4.*dxdxux - 2.*T.^2.*ux.^3.*uy.^4.*dxdxux - 4.*T.^2.*uy.*dxdxuy - 12.*T.^2.*ux.^2.*uy.*dxdxuy - 11.*T.^2.*ux.^4.*uy.*dxdxuy - 3.*T.^2.*ux.^6.*uy.*dxdxuy - 8.*T.^2.*uy.^3.*dxdxuy - 16.*T.^2.*ux.^2.*uy.^3.*dxdxuy - 7.*T.^2.*ux.^4.*uy.^3.*dxdxuy - 5.*T.^2.*uy.^5.*dxdxuy - 5.*T.^2.*ux.^2.*uy.^5.*dxdxuy - T.^2.*uy.^7.*dxdxuy));
    %a = 1;
    
    R2 = (1./(9.*(2 + ux.^2 + uy.^2).^2)).*(8.*pi.^2.*(48.*pi.*T.^2.*ux.*uy.*dyT + 48.*pi.*T.^2.*ux.^3.*uy.*dyT + 12.*pi.*T.^2.*ux.^5.*uy.*dyT + 48.*pi.*T.^2.*ux.*uy.^3.*dyT + 24.*pi.*T.^2.*ux.^3.*uy.^3.*dyT + 12.*pi.*T.^2.*ux.*uy.^5.*dyT + 4.*ux.*dyT.^2 + 4.*ux.^3.*dyT.^2 + ux.^5.*dyT.^2 + 4.*ux.*uy.^2.*dyT.^2 + 2.*ux.^3.*uy.^2.*dyT.^2 + ux.*uy.^4.*dyT.^2 + 16.*pi.*T.^3.*uy.*dyux + 16.*pi.*T.^3.*ux.^2.*uy.*dyux + 4.*pi.*T.^3.*ux.^4.*uy.*dyux + 16.*pi.*T.^3.*uy.^3.*dyux + 8.*pi.*T.^3.*ux.^2.*uy.^3.*dyux + 4.*pi.*T.^3.*uy.^5.*dyux - 4.*T.*dyT.*dyux - 4.*T.*ux.^2.*dyT.*dyux - T.*ux.^4.*dyT.*dyux - 4.*T.*uy.^2.*dyT.*dyux - 2.*T.*ux.^2.*uy.^2.*dyT.*dyux - T.*uy.^4.*dyT.*dyux + 16.*pi.*T.^3.*ux.*dyuy + 16.*pi.*T.^3.*ux.^3.*dyuy + 4.*pi.*T.^3.*ux.^5.*dyuy + 16.*pi.*T.^3.*ux.*uy.^2.*dyuy + 8.*pi.*T.^3.*ux.^3.*uy.^2.*dyuy + 4.*pi.*T.^3.*ux.*uy.^4.*dyuy + 4.*T.*ux.*dydyT + 4.*T.*ux.^3.*dydyT + T.*ux.^5.*dydyT + 4.*T.*ux.*uy.^2.*dydyT + 2.*T.*ux.^3.*uy.^2.*dydyT + T.*ux.*uy.^4.*dydyT - 4.*T.^2.*dydyux - 4.*T.^2.*ux.^2.*dydyux - T.^2.*ux.^4.*dydyux - 4.*T.^2.*uy.^2.*dydyux - 2.*T.^2.*ux.^2.*uy.^2.*dydyux - T.^2.*uy.^4.*dydyux + 16.*pi.*T.^2.*dxT + 64.*pi.*T.^2.*ux.^2.*dxT + 52.*pi.*T.^2.*ux.^4.*dxT + 12.*pi.*T.^2.*ux.^6.*dxT + 16.*pi.*T.^2.*uy.^2.*dxT + 56.*pi.*T.^2.*ux.^2.*uy.^2.*dxT + 24.*pi.*T.^2.*ux.^4.*uy.^2.*dxT + 4.*pi.*T.^2.*uy.^4.*dxT + 12.*pi.*T.^2.*ux.^2.*uy.^4.*dxT + 2.*ux.^2.*uy.*dyT.*dxT + ux.^4.*uy.*dyT.*dxT + 2.*uy.^3.*dyT.*dxT + 2.*ux.^2.*uy.^3.*dyT.*dxT + uy.^5.*dyT.*dxT - 8.*T.*ux.*uy.*dyux.*dxT - 4.*T.*ux.^3.*uy.*dyux.*dxT - 4.*T.*ux.*uy.^3.*dyux.*dxT + 12.*T.*dyuy.*dxT + 16.*T.*ux.^2.*dyuy.*dxT + 5.*T.*ux.^4.*dyuy.*dxT + 8.*T.*uy.^2.*dyuy.*dxT + 6.*T.*ux.^2.*uy.^2.*dyuy.*dxT + T.*uy.^4.*dyuy.*dxT + 4.*ux.*dxT.^2 + 6.*ux.^3.*dxT.^2 + 2.*ux.^5.*dxT.^2 + 6.*ux.*uy.^2.*dxT.^2 + 4.*ux.^3.*uy.^2.*dxT.^2 + 2.*ux.*uy.^4.*dxT.^2 + 32.*pi.*T.^3.*ux.*dxux + 32.*pi.*T.^3.*ux.^3.*dxux + 8.*pi.*T.^3.*ux.^5.*dxux + 32.*pi.*T.^3.*ux.*uy.^2.*dxux + 16.*pi.*T.^3.*ux.^3.*uy.^2.*dxux + 8.*pi.*T.^3.*ux.*uy.^4.*dxux + 4.*T.*ux.*uy.*dyT.*dxux - 4.*T.^2.*uy.*dyux.*dxux + 2.*T.^2.*ux.^2.*uy.*dyux.*dxux - 2.*T.^2.*uy.^3.*dyux.*dxux + 4.*T.^2.*ux.*dyuy.*dxux + 4.*T.^2.*ux.*uy.^2.*dyuy.*dxux - 4.*T.*dxT.*dxux - 2.*T.*ux.^2.*dxT.*dxux - 2.*T.*ux.^4.*dxT.*dxux + 2.*T.*uy.^2.*dxT.*dxux + 2.*T.*uy.^4.*dxT.*dxux - 4.*T.^2.*ux.*dxux.^2 - 4.*T.^2.*ux.*uy.^2.*dxux.^2 - 12.*T.*dyT.*dxuy - 10.*T.*ux.^2.*dyT.*dxuy - 2.*T.*ux.^4.*dyT.*dxuy - 6.*T.*uy.^2.*dyT.*dxuy - 4.*T.*ux.^2.*uy.^2.*dyT.*dxuy - 2.*T.*uy.^4.*dyT.*dxuy - 4.*T.^2.*ux.*dyux.*dxuy - 2.*T.^2.*ux.^3.*dyux.*dxuy + 2.*T.^2.*ux.*uy.^2.*dyux.*dxuy - 4.*T.^2.*uy.*dyuy.*dxuy - 4.*T.^2.*ux.^2.*uy.*dyuy.*dxuy - 4.*T.*ux.*uy.*dxT.*dxuy - 4.*T.*ux.^3.*uy.*dxT.*dxuy - 4.*T.*ux.*uy.^3.*dxT.*dxuy + 6.*T.^2.*ux.^2.*uy.*dxux.*dxuy - 2.*T.^2.*uy.^3.*dxux.*dxuy - 4.*T.^2.*ux.*dxuy.^2 - 2.*T.^2.*ux.^3.*dxuy.^2 + 2.*T.^2.*ux.*uy.^2.*dxuy.^2 + 2.*T.*ux.^2.*uy.*dxdyT + T.*ux.^4.*uy.*dxdyT + 2.*T.*uy.^3.*dxdyT + 2.*T.*ux.^2.*uy.^3.*dxdyT + T.*uy.^5.*dxdyT - 4.*T.^2.*ux.*uy.*dxdyux - 2.*T.^2.*ux.^3.*uy.*dxdyux - 2.*T.^2.*ux.*uy.^3.*dxdyux + 2.*T.^2.*ux.^2.*dxdyuy + T.^2.*ux.^4.*dxdyuy - 2.*T.^2.*uy.^2.*dxdyuy - T.^2.*uy.^4.*dxdyuy + 4.*T.*ux.*dxdxT + 6.*T.*ux.^3.*dxdxT + 2.*T.*ux.^5.*dxdxT + 6.*T.*ux.*uy.^2.*dxdxT + 4.*T.*ux.^3.*uy.^2.*dxdxT + 2.*T.*ux.*uy.^4.*dxdxT - 4.*T.^2.*dxdxux - 6.*T.^2.*ux.^2.*dxdxux - 2.*T.^2.*ux.^4.*dxdxux - 2.*T.^2.*uy.^2.*dxdxux - 2.*T.^2.*ux.^2.*uy.^2.*dxdxux - 4.*T.^2.*ux.*uy.*dxdxuy - 2.*T.^2.*ux.^3.*uy.*dxdxuy - 2.*T.^2.*ux.*uy.^3.*dxdxuy)); 
    %R2 = R2 + a * (    dxdx(ux.*dxux + uy.*dyux) + dydy(ux.*dxux + uy.*dyux)  ) + a * (    dxdx(ux.*dxux + uy.*dyux) + dx(dy(ux.*dxuy + uy.*dyuy))  );
    %R2 = R2 + a * ( ux.*dx(dxux + dyuy) + uy.*dy(dxux + dyuy));

    R3 = (1./(9.*(2 + ux.^2 + uy.^2).^2)).*(8.*pi.^2.*(16.*pi.*T.^2.*dyT + 16.*pi.*T.^2.*ux.^2.*dyT + 4.*pi.*T.^2.*ux.^4.*dyT + 64.*pi.*T.^2.*uy.^2.*dyT + 56.*pi.*T.^2.*ux.^2.*uy.^2.*dyT + 12.*pi.*T.^2.*ux.^4.*uy.^2.*dyT + 52.*pi.*T.^2.*uy.^4.*dyT + 24.*pi.*T.^2.*ux.^2.*uy.^4.*dyT + 12.*pi.*T.^2.*uy.^6.*dyT + 4.*uy.*dyT.^2 + 6.*ux.^2.*uy.*dyT.^2 + 2.*ux.^4.*uy.*dyT.^2 + 6.*uy.^3.*dyT.^2 + 4.*ux.^2.*uy.^3.*dyT.^2 + 2.*uy.^5.*dyT.^2 - 4.*T.*ux.*uy.*dyT.*dyux - 4.*T.*ux.^3.*uy.*dyT.*dyux - 4.*T.*ux.*uy.^3.*dyT.*dyux - 4.*T.^2.*uy.*dyux.^2 + 2.*T.^2.*ux.^2.*uy.*dyux.^2 - 2.*T.^2.*uy.^3.*dyux.^2 + 32.*pi.*T.^3.*uy.*dyuy + 32.*pi.*T.^3.*ux.^2.*uy.*dyuy + 8.*pi.*T.^3.*ux.^4.*uy.*dyuy + 32.*pi.*T.^3.*uy.^3.*dyuy + 16.*pi.*T.^3.*ux.^2.*uy.^3.*dyuy + 8.*pi.*T.^3.*uy.^5.*dyuy - 4.*T.*dyT.*dyuy + 2.*T.*ux.^2.*dyT.*dyuy + 2.*T.*ux.^4.*dyT.*dyuy - 2.*T.*uy.^2.*dyT.*dyuy - 2.*T.*uy.^4.*dyT.*dyuy - 2.*T.^2.*ux.^3.*dyux.*dyuy + 6.*T.^2.*ux.*uy.^2.*dyux.*dyuy - 4.*T.^2.*uy.*dyuy.^2 - 4.*T.^2.*ux.^2.*uy.*dyuy.^2 + 4.*T.*uy.*dydyT + 6.*T.*ux.^2.*uy.*dydyT + 2.*T.*ux.^4.*uy.*dydyT + 6.*T.*uy.^3.*dydyT + 4.*T.*ux.^2.*uy.^3.*dydyT + 2.*T.*uy.^5.*dydyT - 4.*T.^2.*ux.*uy.*dydyux - 2.*T.^2.*ux.^3.*uy.*dydyux - 2.*T.^2.*ux.*uy.^3.*dydyux - 4.*T.^2.*dydyuy - 2.*T.^2.*ux.^2.*dydyuy - 6.*T.^2.*uy.^2.*dydyuy - 2.*T.^2.*ux.^2.*uy.^2.*dydyuy - 2.*T.^2.*uy.^4.*dydyuy + 48.*pi.*T.^2.*ux.*uy.*dxT + 48.*pi.*T.^2.*ux.^3.*uy.*dxT + 12.*pi.*T.^2.*ux.^5.*uy.*dxT + 48.*pi.*T.^2.*ux.*uy.^3.*dxT + 24.*pi.*T.^2.*ux.^3.*uy.^3.*dxT + 12.*pi.*T.^2.*ux.*uy.^5.*dxT + 2.*ux.^3.*dyT.*dxT + ux.^5.*dyT.*dxT + 2.*ux.*uy.^2.*dyT.*dxT + 2.*ux.^3.*uy.^2.*dyT.*dxT + ux.*uy.^4.*dyT.*dxT - 12.*T.*dyux.*dxT - 6.*T.*ux.^2.*dyux.*dxT - 2.*T.*ux.^4.*dyux.*dxT - 10.*T.*uy.^2.*dyux.*dxT - 4.*T.*ux.^2.*uy.^2.*dyux.*dxT - 2.*T.*uy.^4.*dyux.*dxT + 4.*T.*ux.*uy.*dyuy.*dxT + 4.*uy.*dxT.^2 + 4.*ux.^2.*uy.*dxT.^2 + ux.^4.*uy.*dxT.^2 + 4.*uy.^3.*dxT.^2 + 2.*ux.^2.*uy.^3.*dxT.^2 + uy.^5.*dxT.^2 + 16.*pi.*T.^3.*uy.*dxux + 16.*pi.*T.^3.*ux.^2.*uy.*dxux + 4.*pi.*T.^3.*ux.^4.*uy.*dxux + 16.*pi.*T.^3.*uy.^3.*dxux + 8.*pi.*T.^3.*ux.^2.*uy.^3.*dxux + 4.*pi.*T.^3.*uy.^5.*dxux + 12.*T.*dyT.*dxux + 8.*T.*ux.^2.*dyT.*dxux + T.*ux.^4.*dyT.*dxux + 16.*T.*uy.^2.*dyT.*dxux + 6.*T.*ux.^2.*uy.^2.*dyT.*dxux + 5.*T.*uy.^4.*dyT.*dxux - 4.*T.^2.*ux.*dyux.*dxux - 4.*T.^2.*ux.*uy.^2.*dyux.*dxux + 4.*T.^2.*uy.*dyuy.*dxux + 4.*T.^2.*ux.^2.*uy.*dyuy.*dxux + 16.*pi.*T.^3.*ux.*dxuy + 16.*pi.*T.^3.*ux.^3.*dxuy + 4.*pi.*T.^3.*ux.^5.*dxuy + 16.*pi.*T.^3.*ux.*uy.^2.*dxuy + 8.*pi.*T.^3.*ux.^3.*uy.^2.*dxuy + 4.*pi.*T.^3.*ux.*uy.^4.*dxuy - 8.*T.*ux.*uy.*dyT.*dxuy - 4.*T.*ux.^3.*uy.*dyT.*dxuy - 4.*T.*ux.*uy.^3.*dyT.*dxuy - 4.*T.^2.*uy.*dyux.*dxuy + 2.*T.^2.*ux.^2.*uy.*dyux.*dxuy - 2.*T.^2.*uy.^3.*dyux.*dxuy - 4.*T.^2.*ux.*dyuy.*dxuy - 2.*T.^2.*ux.^3.*dyuy.*dxuy + 2.*T.^2.*ux.*uy.^2.*dyuy.*dxuy - 4.*T.*dxT.*dxuy - 4.*T.*ux.^2.*dxT.*dxuy - T.*ux.^4.*dxT.*dxuy - 4.*T.*uy.^2.*dxT.*dxuy - 2.*T.*ux.^2.*uy.^2.*dxT.*dxuy - T.*uy.^4.*dxT.*dxuy + 2.*T.*ux.^3.*dxdyT + T.*ux.^5.*dxdyT + 2.*T.*ux.*uy.^2.*dxdyT + 2.*T.*ux.^3.*uy.^2.*dxdyT + T.*ux.*uy.^4.*dxdyT - 2.*T.^2.*ux.^2.*dxdyux - T.^2.*ux.^4.*dxdyux + 2.*T.^2.*uy.^2.*dxdyux + T.^2.*uy.^4.*dxdyux - 4.*T.^2.*ux.*uy.*dxdyuy - 2.*T.^2.*ux.^3.*uy.*dxdyuy - 2.*T.^2.*ux.*uy.^3.*dxdyuy + 4.*T.*uy.*dxdxT + 4.*T.*ux.^2.*uy.*dxdxT + T.*ux.^4.*uy.*dxdxT + 4.*T.*uy.^3.*dxdxT + 2.*T.*ux.^2.*uy.^3.*dxdxT + T.*uy.^5.*dxdxT - 4.*T.^2.*dxdxuy - 4.*T.^2.*ux.^2.*dxdxuy - T.^2.*ux.^4.*dxdxuy - 4.*T.^2.*uy.^2.*dxdxuy - 2.*T.^2.*ux.^2.*uy.^2.*dxdxuy - T.^2.*uy.^4.*dxdxuy));
    %R3 = R3 + a * (    dx(dy(ux.*dxux + uy.*dyux)) + dydy(ux.*dxuy + uy.*dyuy)  ) + a * (    dxdx(ux.*dxuy + uy.*dyuy) + dydy(ux.*dxuy + uy.*dyuy)  );
    %R2 = R2 + a * ( ux.*dx(dxux + dyuy) + uy.*dy(dxux + dyuy));

    R = [R1 R2 R3];


    [dt, fval, exitflag, output, jacobian] = fsolve(@eq, zeros(GRID_SIZE, GRID_SIZE*3), optimset('Jacobian', 'on', 'DerivativeCheck', 'off', 'Display', 'off', 'Diagnostics', 'off'));

    dtT  = dt(:,            1:  GRID_SIZE);
    dtux = dt(:,  GRID_SIZE+1:2*GRID_SIZE);
    dtuy = dt(:,2*GRID_SIZE+1:3*GRID_SIZE);

  



    ux = ux + DT * dtux;
    uy = uy + DT * dtuy;
    T  = T  + DT * dtT;


    dxux = dx(ux);
    dyux = dy(ux);
    dxuy = dx(uy);
    dyuy = dy(uy);
    dxT = dx(T);
    dyT = dy(T);

    dxdxux = dxdx(ux);
    dxdyux = dx(dy(ux));
    dydyux = dydy(ux);
    dxdxuy = dxdx(uy);
    dxdyuy = dx(dy(uy));
    dydyuy = dydy(uy);
    dxdxT = dxdx(T);
    dxdyT = dx(dy(T));
    dydyT = dydy(T);



    phydro();
    pause(0.003);
  end

















    phydro();
    %pause(0.033);

  function [res J] = eq(g)
    dtT  = g(:,            1:  GRID_SIZE);
    dtux = g(:,  GRID_SIZE+1:2*GRID_SIZE);
    dtuy = g(:,2*GRID_SIZE+1:3*GRID_SIZE);

    dxdtT = dx(dtT);
    dydtT = dy(dtT);
    dxdtux = dx(dtux);
    dydtux = dy(dtux);
    dxdtuy = dx(dtuy);
    dydtuy = dy(dtuy);

    L1 = (1./(9.*(2 + ux.^2 + uy.^2).^2)).*(16.*pi.^2.*(16.*pi.*T.^2.*dtT + 40.*pi.*T.^2.*ux.^2.*dtT + 28.*pi.*T.^2.*ux.^4.*dtT + 6.*pi.*T.^2.*ux.^6.*dtT + 40.*pi.*T.^2.*uy.^2.*dtT + 56.*pi.*T.^2.*ux.^2.*uy.^2.*dtT + 18.*pi.*T.^2.*ux.^4.*uy.^2.*dtT + 28.*pi.*T.^2.*uy.^4.*dtT + 18.*pi.*T.^2.*ux.^2.*uy.^4.*dtT + 6.*pi.*T.^2.*uy.^6.*dtT + 2.*ux.^2.*uy.*dyT.*dtT + ux.^4.*uy.*dyT.*dtT + 2.*uy.^3.*dyT.*dtT + 2.*ux.^2.*uy.^3.*dyT.*dtT + uy.^5.*dyT.*dtT - 8.*T.*ux.*uy.*dyux.*dtT - 4.*T.*ux.^3.*uy.*dyux.*dtT - 4.*T.*ux.*uy.^3.*dyux.*dtT + 4.*T.*ux.^2.*dyuy.*dtT + 2.*T.*ux.^4.*dyuy.*dtT - 4.*T.*uy.^2.*dyuy.*dtT - 2.*T.*uy.^4.*dyuy.*dtT + 2.*ux.^3.*dxT.*dtT + ux.^5.*dxT.*dtT + 2.*ux.*uy.^2.*dxT.*dtT + 2.*ux.^3.*uy.^2.*dxT.*dtT + ux.*uy.^4.*dxT.*dtT - 4.*T.*ux.^2.*dxux.*dtT - 2.*T.*ux.^4.*dxux.*dtT + 4.*T.*uy.^2.*dxux.*dtT + 2.*T.*uy.^4.*dxux.*dtT - 8.*T.*ux.*uy.*dxuy.*dtT - 4.*T.*ux.^3.*uy.*dxuy.*dtT - 4.*T.*ux.*uy.^3.*dxuy.*dtT + 16.*pi.*T.^3.*ux.*dtux + 16.*pi.*T.^3.*ux.^3.*dtux + 4.*pi.*T.^3.*ux.^5.*dtux + 16.*pi.*T.^3.*ux.*uy.^2.*dtux + 8.*pi.*T.^3.*ux.^3.*uy.^2.*dtux + 4.*pi.*T.^3.*ux.*uy.^4.*dtux + 4.*T.*ux.*uy.*dyT.*dtux - 4.*T.^2.*uy.*dyux.*dtux + 2.*T.^2.*ux.^2.*uy.*dyux.*dtux - 2.*T.^2.*uy.^3.*dyux.*dtux + 4.*T.^2.*ux.*dyuy.*dtux + 4.*T.^2.*ux.*uy.^2.*dyuy.*dtux + 6.*T.*ux.^2.*dxT.*dtux + T.*ux.^4.*dxT.*dtux + 2.*T.*uy.^2.*dxT.*dtux + 2.*T.*ux.^2.*uy.^2.*dxT.*dtux + T.*uy.^4.*dxT.*dtux - 4.*T.^2.*ux.*dxux.*dtux - 4.*T.^2.*ux.*uy.^2.*dxux.*dtux - 4.*T.^2.*uy.*dxuy.*dtux + 2.*T.^2.*ux.^2.*uy.*dxuy.*dtux - 2.*T.^2.*uy.^3.*dxuy.*dtux + 16.*pi.*T.^3.*uy.*dtuy + 16.*pi.*T.^3.*ux.^2.*uy.*dtuy + 4.*pi.*T.^3.*ux.^4.*uy.*dtuy + 16.*pi.*T.^3.*uy.^3.*dtuy + 8.*pi.*T.^3.*ux.^2.*uy.^3.*dtuy + 4.*pi.*T.^3.*uy.^5.*dtuy + 2.*T.*ux.^2.*dyT.*dtuy + T.*ux.^4.*dyT.*dtuy + 6.*T.*uy.^2.*dyT.*dtuy + 2.*T.*ux.^2.*uy.^2.*dyT.*dtuy + T.*uy.^4.*dyT.*dtuy - 4.*T.^2.*ux.*dyux.*dtuy - 2.*T.^2.*ux.^3.*dyux.*dtuy + 2.*T.^2.*ux.*uy.^2.*dyux.*dtuy - 4.*T.^2.*uy.*dyuy.*dtuy - 4.*T.^2.*ux.^2.*uy.*dyuy.*dtuy + 4.*T.*ux.*uy.*dxT.*dtuy + 4.*T.^2.*uy.*dxux.*dtuy + 4.*T.^2.*ux.^2.*uy.*dxux.*dtuy - 4.*T.^2.*ux.*dxuy.*dtuy - 2.*T.^2.*ux.^3.*dxuy.*dtuy + 2.*T.^2.*ux.*uy.^2.*dxuy.*dtuy + 2.*T.*ux.^2.*uy.*dydtT + T.*ux.^4.*uy.*dydtT + 2.*T.*uy.^3.*dydtT + 2.*T.*ux.^2.*uy.^3.*dydtT + T.*uy.^5.*dydtT - 4.*T.^2.*ux.*uy.*dydtux - 2.*T.^2.*ux.^3.*uy.*dydtux - 2.*T.^2.*ux.*uy.^3.*dydtux + 2.*T.^2.*ux.^2.*dydtuy + T.^2.*ux.^4.*dydtuy - 2.*T.^2.*uy.^2.*dydtuy - T.^2.*uy.^4.*dydtuy + 2.*T.*ux.^3.*dxdtT + T.*ux.^5.*dxdtT + 2.*T.*ux.*uy.^2.*dxdtT + 2.*T.*ux.^3.*uy.^2.*dxdtT + T.*ux.*uy.^4.*dxdtT - 2.*T.^2.*ux.^2.*dxdtux - T.^2.*ux.^4.*dxdtux + 2.*T.^2.*uy.^2.*dxdtux + T.^2.*uy.^4.*dxdtux - 4.*T.^2.*ux.*uy.*dxdtuy - 2.*T.^2.*ux.^3.*uy.*dxdtuy - 2.*T.^2.*ux.*uy.^3.*dxdtuy));

    L2 = (1./(9.*(1 + ux.^2 + uy.^2).^(3./2).*(2 + ux.^2 + uy.^2).^2)).*(8.* pi.^2.*(48.*pi.*T.^2.*ux.* dtT + 144.*pi.*T.^2.*ux.^3.* dtT + 156.*pi.*T.^2.*ux.^5.* dtT + 72.*pi.*T.^2.*ux.^7.* dtT + 12.*pi.*T.^2.*ux.^9.* dtT + 144.*pi.*T.^2.*ux.*uy.^2.* dtT + 312.*pi.*T.^2.*ux.^3.*uy.^2.* dtT + 216.*pi.*T.^2.*ux.^5.*uy.^2.* dtT + 48.*pi.*T.^2.*ux.^7.*uy.^2.* dtT + 156.*pi.*T.^2.*ux.*uy.^4.* dtT + 216.*pi.*T.^2.*ux.^3.*uy.^4.* dtT + 72.*pi.*T.^2.*ux.^5.*uy.^4.* dtT + 72.*pi.*T.^2.*ux.*uy.^6.* dtT + 48.*pi.*T.^2.*ux.^3.*uy.^6.* dtT + 12.*pi.*T.^2.*ux.*uy.^8.* dtT + 2.*ux.^3.*uy.*dyT.* dtT + 3.*ux.^5.*uy.*dyT.* dtT + ux.^7.*uy.*dyT.* dtT + 2.*ux.*uy.^3.*dyT.* dtT + 6.*ux.^3.*uy.^3.*dyT.* dtT + 3.*ux.^5.*uy.^3.*dyT.* dtT + 3.*ux.*uy.^5.*dyT.* dtT + 3.*ux.^3.*uy.^5.*dyT.* dtT + ux.*uy.^7.*dyT.* dtT - 8.*T.*uy.*dyux.* dtT - 24.*T.*ux.^2.*uy.* dyux.* dtT - 22.*T.*ux.^4.*uy.* dyux.* dtT - 6.*T.*ux.^6.*uy.* dyux.* dtT - 16.*T.*uy.^3.*dyux.* dtT - 32.*T.*ux.^2.*uy.^3.* dyux.* dtT - 14.*T.*ux.^4.*uy.^3.* dyux.* dtT - 10.*T.*uy.^5.*dyux.* dtT - 10.*T.*ux.^2.*uy.^5.* dyux.* dtT - 2.*T.*uy.^7.*dyux.* dtT + 8.*T.*ux.*dyuy.* dtT + 20.*T.*ux.^3.*dyuy.* dtT + 16.*T.*ux.^5.*dyuy.* dtT + 4.*T.*ux.^7.*dyuy.* dtT + 12.*T.*ux.*uy.^2.* dyuy.* dtT + 20.*T.*ux.^3.*uy.^2.* dyuy.* dtT + 8.*T.*ux.^5.*uy.^2.* dyuy.* dtT + 4.*T.*ux.*uy.^4.* dyuy.* dtT + 4.*T.*ux.^3.*uy.^4.* dyuy.* dtT + 4.*ux.^2.*dxT.* dtT + 10.*ux.^4.*dxT.* dtT + 8.*ux.^6.*dxT.* dtT + 2.*ux.^8.*dxT.* dtT + 4.*uy.^2.*dxT.* dtT + 18.*ux.^2.*uy.^2.*dxT.* dtT + 21.*ux.^4.*uy.^2.*dxT.* dtT + 7.*ux.^6.*uy.^2.*dxT.* dtT + 8.*uy.^4.*dxT.* dtT + 18.*ux.^2.*uy.^4.*dxT.* dtT + 9.*ux.^4.*uy.^4.*dxT.* dtT + 5.*uy.^6.*dxT.* dtT + 5.*ux.^2.*uy.^6.*dxT.* dtT + uy.^8.*dxT.* dtT - 8.*T.*ux.*dxux.* dtT - 20.*T.*ux.^3.*dxux.* dtT - 16.*T.*ux.^5.*dxux.* dtT - 4.*T.*ux.^7.*dxux.* dtT - 12.*T.*ux.*uy.^2.* dxux.* dtT - 20.*T.*ux.^3.*uy.^2.* dxux.* dtT - 8.*T.*ux.^5.*uy.^2.* dxux.* dtT - 4.*T.*ux.*uy.^4.* dxux.* dtT - 4.*T.*ux.^3.*uy.^4.* dxux.* dtT - 8.*T.*uy.*dxuy.* dtT - 24.*T.*ux.^2.*uy.* dxuy.* dtT - 22.*T.*ux.^4.*uy.* dxuy.* dtT - 6.*T.*ux.^6.*uy.* dxuy.* dtT - 16.*T.*uy.^3.*dxuy.* dtT - 32.*T.*ux.^2.*uy.^3.* dxuy.* dtT - 14.*T.*ux.^4.*uy.^3.* dxuy.* dtT - 10.*T.*uy.^5.*dxuy.* dtT - 10.*T.*ux.^2.*uy.^5.* dxuy.* dtT - 2.*T.*uy.^7.*dxuy.* dtT + 16.*pi.*T.^3.*dtux + 64.*pi.*T.^3.*ux.^2.* dtux + 84.*pi.*T.^3.*ux.^4.* dtux + 44.*pi.*T.^3.*ux.^6.* dtux + 8.*pi.*T.^3.*ux.^8.* dtux + 48.*pi.*T.^3.*uy.^2.* dtux + 136.*pi.*T.^3.*ux.^2.*uy.^2.* dtux + 112.*pi.*T.^3.*ux.^4.*uy.^2.* dtux + 28.*pi.*T.^3.*ux.^6.*uy.^2.* dtux + 52.*pi.*T.^3.*uy.^4.* dtux + 92.*pi.*T.^3.*ux.^2.*uy.^4.* dtux + 36.*pi.*T.^3.*ux.^4.*uy.^4.* dtux + 24.*pi.*T.^3.*uy.^6.* dtux + 20.*pi.*T.^3.*ux.^2.*uy.^6.* dtux + 4.*pi.*T.^3.*uy.^8.* dtux + 6.*T.*ux.^2.*uy.* dyT.* dtux + 5.*T.*ux.^4.*uy.* dyT.* dtux + 2.*T.*uy.^3.*dyT.* dtux + 8.*T.*ux.^2.*uy.^3.* dyT.* dtux + T.*ux.^4.*uy.^3.* dyT.* dtux + 3.*T.*uy.^5.*dyT.* dtux + 2.*T.*ux.^2.*uy.^5.* dyT.* dtux + T.*uy.^7.*dyT.* dtux - 4.*T.^2.*ux.*uy.* dyux.* dtux + 3.*T.^2.*ux.^5.*uy.* dyux.* dtux - 8.*T.^2.*ux.*uy.^3.* dyux.* dtux - 3.*T.^2.*ux.*uy.^5.* dyux.* dtux + 4.*T.^2.*dyuy.* dtux + 10.*T.^2.*ux.^2.*dyuy.* dtux + 6.*T.^2.*ux.^4.*dyuy.* dtux + 6.*T.^2.*uy.^2.*dyuy.* dtux + 16.*T.^2.*ux.^2.*uy.^2.* dyuy.* dtux + 6.*T.^2.*ux.^4.*uy.^2.* dyuy.* dtux + 2.*T.^2.*uy.^4.*dyuy.* dtux + 6.*T.^2.*ux.^2.*uy.^4.* dyuy.* dtux + 8.*T.*ux.*dxT.* dtux + 20.*T.*ux.^3.*dxT.* dtux + 14.*T.*ux.^5.*dxT.* dtux + 2.*T.*ux.^7.*dxT.* dtux + 16.*T.*ux.*uy.^2.* dxT.* dtux + 26.*T.*ux.^3.*uy.^2.* dxT.* dtux + 7.*T.*ux.^5.*uy.^2.* dxT.* dtux + 12.*T.*ux.*uy.^4.* dxT.* dtux + 8.*T.*ux.^3.*uy.^4.* dxT.* dtux + 3.*T.*ux.*uy.^6.* dxT.* dtux - 4.*T.^2.*dxux.* dtux - 10.*T.^2.*ux.^2.*dxux.* dtux - 6.*T.^2.*ux.^4.*dxux.* dtux - 6.*T.^2.*uy.^2.*dxux.* dtux - 16.*T.^2.*ux.^2.*uy.^2.* dxux.* dtux - 6.*T.^2.*ux.^4.*uy.^2.* dxux.* dtux - 2.*T.^2.*uy.^4.*dxux.* dtux - 6.*T.^2.*ux.^2.*uy.^4.* dxux.* dtux - 4.*T.^2.*ux.*uy.* dxuy.* dtux + 3.*T.^2.*ux.^5.*uy.* dxuy.* dtux - 8.*T.^2.*ux.*uy.^3.* dxuy.* dtux - 3.*T.^2.*ux.*uy.^5.* dxuy.* dtux + 16.*pi.*T.^3.*ux.*uy.* dtuy + 32.*pi.*T.^3.*ux.^3.*uy.* dtuy + 20.*pi.*T.^3.*ux.^5.*uy.* dtuy + 4.*pi.*T.^3.*ux.^7.*uy.* dtuy + 32.*pi.*T.^3.*ux.*uy.^3.* dtuy + 40.*pi.*T.^3.*ux.^3.*uy.^3.* dtuy + 12.*pi.*T.^3.*ux.^5.*uy.^3.* dtuy + 20.*pi.*T.^3.*ux.*uy.^5.* dtuy + 12.*pi.*T.^3.*ux.^3.*uy.^5.* dtuy + 4.*pi.*T.^3.*ux.*uy.^7.* dtuy + 2.*T.*ux.^3.*dyT.* dtuy + 3.*T.*ux.^5.*dyT.* dtuy + T.*ux.^7.*dyT.* dtuy + 6.*T.*ux.*uy.^2.* dyT.* dtuy + 8.*T.*ux.^3.*uy.^2.* dyT.* dtuy + 2.*T.*ux.^5.*uy.^2.* dyT.* dtuy + 5.*T.*ux.*uy.^4.* dyT.* dtuy + T.*ux.^3.*uy.^4.* dyT.* dtuy - 4.*T.^2.*dyux.* dtuy - 12.*T.^2.*ux.^2.*dyux.* dtuy - 11.*T.^2.*ux.^4.*dyux.* dtuy - 3.*T.^2.*ux.^6.*dyux.* dtuy - 4.*T.^2.*uy.^2.*dyux.* dtuy - 4.*T.^2.*ux.^2.*uy.^2.* dyux.* dtuy - T.^2.*uy.^4.*dyux.* dtuy + 3.*T.^2.*ux.^2.*uy.^4.* dyux.* dtuy - 8.*T.^2.*ux.*uy.* dyuy.* dtuy - 14.*T.^2.*ux.^3.*uy.* dyuy.* dtuy - 6.*T.^2.*ux.^5.*uy.* dyuy.* dtuy - 6.*T.^2.*ux.*uy.^3.* dyuy.* dtuy - 6.*T.^2.*ux.^3.*uy.^3.* dyuy.* dtuy + 8.*T.*uy.*dxT.* dtuy + 16.*T.*ux.^2.*uy.* dxT.* dtuy + 8.*T.*ux.^4.*uy.* dxT.* dtuy + 12.*T.*uy.^3.*dxT.* dtuy + 14.*T.*ux.^2.*uy.^3.* dxT.* dtuy + T.*ux.^4.*uy.^3.* dxT.* dtuy + 6.*T.*uy.^5.*dxT.* dtuy + 2.*T.*ux.^2.*uy.^5.* dxT.* dtuy + T.*uy.^7.*dxT.* dtuy + 8.*T.^2.*ux.*uy.* dxux.* dtuy + 14.*T.^2.*ux.^3.*uy.* dxux.* dtuy + 6.*T.^2.*ux.^5.*uy.* dxux.* dtuy + 6.*T.^2.*ux.*uy.^3.* dxux.* dtuy + 6.*T.^2.*ux.^3.*uy.^3.* dxux.* dtuy - 4.*T.^2.*dxuy.* dtuy - 12.*T.^2.*ux.^2.*dxuy.* dtuy - 11.*T.^2.*ux.^4.*dxuy.* dtuy - 3.*T.^2.*ux.^6.*dxuy.* dtuy - 4.*T.^2.*uy.^2.*dxuy.* dtuy - 4.*T.^2.*ux.^2.*uy.^2.* dxuy.* dtuy - T.^2.*uy.^4.*dxuy.* dtuy + 3.*T.^2.*ux.^2.*uy.^4.* dxuy.* dtuy + 2.*T.*ux.^3.*uy.* dydtT + 3.*T.*ux.^5.*uy.* dydtT + T.*ux.^7.*uy.* dydtT + 2.*T.*ux.*uy.^3.* dydtT + 6.*T.*ux.^3.*uy.^3.* dydtT + 3.*T.*ux.^5.*uy.^3.* dydtT + 3.*T.*ux.*uy.^5.* dydtT + 3.*T.*ux.^3.*uy.^5.* dydtT + T.*ux.*uy.^7.* dydtT - 4.*T.^2.*uy.*dydtux - 12.*T.^2.*ux.^2.*uy.* dydtux - 11.*T.^2.*ux.^4.*uy.* dydtux - 3.*T.^2.*ux.^6.*uy.* dydtux - 8.*T.^2.*uy.^3.*dydtux - 16.*T.^2.*ux.^2.*uy.^3.* dydtux - 7.*T.^2.*ux.^4.*uy.^3.* dydtux - 5.*T.^2.*uy.^5.*dydtux - 5.*T.^2.*ux.^2.*uy.^5.* dydtux - T.^2.*uy.^7.*dydtux + 4.*T.^2.*ux.*dydtuy + 10.*T.^2.*ux.^3.* dydtuy + 8.*T.^2.*ux.^5.*dydtuy + 2.*T.^2.*ux.^7.* dydtuy + 6.*T.^2.*ux.*uy.^2.* dydtuy + 10.*T.^2.*ux.^3.*uy.^2.* dydtuy + 4.*T.^2.*ux.^5.*uy.^2.* dydtuy + 2.*T.^2.*ux.*uy.^4.* dydtuy + 2.*T.^2.*ux.^3.*uy.^4.* dydtuy + 4.*T.*ux.^2.*dxdtT + 10.*T.*ux.^4.*dxdtT + 8.*T.*ux.^6.*dxdtT + 2.*T.*ux.^8.*dxdtT + 4.*T.*uy.^2.*dxdtT + 18.*T.*ux.^2.*uy.^2.* dxdtT + 21.*T.*ux.^4.*uy.^2.* dxdtT + 7.*T.*ux.^6.*uy.^2.* dxdtT + 8.*T.*uy.^4.*dxdtT + 18.*T.*ux.^2.*uy.^4.* dxdtT + 9.*T.*ux.^4.*uy.^4.* dxdtT + 5.*T.*uy.^6.*dxdtT + 5.*T.*ux.^2.*uy.^6.* dxdtT + T.*uy.^8.*dxdtT - 4.*T.^2.*ux.*dxdtux - 10.*T.^2.*ux.^3.* dxdtux - 8.*T.^2.*ux.^5.*dxdtux - 2.*T.^2.*ux.^7.*dxdtux - 6.*T.^2.*ux.*uy.^2.* dxdtux - 10.*T.^2.*ux.^3.*uy.^2.* dxdtux - 4.*T.^2.*ux.^5.*uy.^2.* dxdtux - 2.*T.^2.*ux.*uy.^4.* dxdtux - 2.*T.^2.*ux.^3.*uy.^4.* dxdtux - 4.*T.^2.*uy.*dxdtuy - 12.*T.^2.*ux.^2.*uy.* dxdtuy - 11.*T.^2.*ux.^4.*uy.* dxdtuy - 3.*T.^2.*ux.^6.*uy.* dxdtuy - 8.*T.^2.*uy.^3.*dxdtuy - 16.*T.^2.*ux.^2.*uy.^3.* dxdtuy - 7.*T.^2.*ux.^4.*uy.^3.* dxdtuy - 5.*T.^2.*uy.^5.*dxdtuy - 5.*T.^2.*ux.^2.*uy.^5.* dxdtuy - T.^2.*uy.^7.* dxdtuy));

    L3 = (1./(9.*(1 + ux.^2 + uy.^2).^(3./2).*(2 + ux.^2 + uy.^2).^2)).*(8.*pi.^2.*(48.*pi.*T.^2.*uy.*dtT + 144.*pi.*T.^2.*ux.^2.*uy.*dtT + 156.*pi.*T.^2.*ux.^4.*uy.*dtT + 72.*pi.*T.^2.*ux.^6.*uy.*dtT + 12.*pi.*T.^2.*ux.^8.*uy.*dtT + 144.*pi.*T.^2.*uy.^3.*dtT + 312.*pi.*T.^2.*ux.^2.*uy.^3.*dtT + 216.*pi.*T.^2.*ux.^4.*uy.^3.*dtT + 48.*pi.*T.^2.*ux.^6.*uy.^3.*dtT + 156.*pi.*T.^2.*uy.^5.*dtT + 216.*pi.*T.^2.*ux.^2.*uy.^5.*dtT + 72.*pi.*T.^2.*ux.^4.*uy.^5.*dtT + 72.*pi.*T.^2.*uy.^7.*dtT + 48.*pi.*T.^2.*ux.^2.*uy.^7.*dtT + 12.*pi.*T.^2.*uy.^9.*dtT + 4.*ux.^2.*dyT.*dtT + 8.*ux.^4.*dyT.*dtT + 5.*ux.^6.*dyT.*dtT + ux.^8.*dyT.*dtT + 4.*uy.^2.*dyT.*dtT + 18.*ux.^2.*uy.^2.*dyT.*dtT + 18.*ux.^4.*uy.^2.*dyT.*dtT + 5.*ux.^6.*uy.^2.*dyT.*dtT + 10.*uy.^4.*dyT.*dtT + 21.*ux.^2.*uy.^4.*dyT.*dtT + 9.*ux.^4.*uy.^4.*dyT.*dtT + 8.*uy.^6.*dyT.*dtT + 7.*ux.^2.*uy.^6.*dyT.*dtT + 2.*uy.^8.*dyT.*dtT - 8.*T.*ux.*dyux.*dtT - 16.*T.*ux.^3.*dyux.*dtT - 10.*T.*ux.^5.*dyux.*dtT - 2.*T.*ux.^7.*dyux.*dtT - 24.*T.*ux.*uy.^2.*dyux.*dtT - 32.*T.*ux.^3.*uy.^2.*dyux.*dtT - 10.*T.*ux.^5.*uy.^2.*dyux.*dtT - 22.*T.*ux.*uy.^4.*dyux.*dtT - 14.*T.*ux.^3.*uy.^4.*dyux.*dtT - 6.*T.*ux.*uy.^6.*dyux.*dtT - 8.*T.*uy.*dyuy.*dtT - 12.*T.*ux.^2.*uy.*dyuy.*dtT - 4.*T.*ux.^4.*uy.*dyuy.*dtT - 20.*T.*uy.^3.*dyuy.*dtT - 20.*T.*ux.^2.*uy.^3.*dyuy.*dtT - 4.*T.*ux.^4.*uy.^3.*dyuy.*dtT - 16.*T.*uy.^5.*dyuy.*dtT - 8.*T.*ux.^2.*uy.^5.*dyuy.*dtT - 4.*T.*uy.^7.*dyuy.*dtT + 2.*ux.^3.*uy.*dxT.*dtT + 3.*ux.^5.*uy.*dxT.*dtT + ux.^7.*uy.*dxT.*dtT + 2.*ux.*uy.^3.*dxT.*dtT + 6.*ux.^3.*uy.^3.*dxT.*dtT + 3.*ux.^5.*uy.^3.*dxT.*dtT + 3.*ux.*uy.^5.*dxT.*dtT + 3.*ux.^3.*uy.^5.*dxT.*dtT + ux.*uy.^7.*dxT.*dtT + 8.*T.*uy.*dxux.*dtT + 12.*T.*ux.^2.*uy.*dxux.*dtT + 4.*T.*ux.^4.*uy.*dxux.*dtT + 20.*T.*uy.^3.*dxux.*dtT + 20.*T.*ux.^2.*uy.^3.*dxux.*dtT + 4.*T.*ux.^4.*uy.^3.*dxux.*dtT + 16.*T.*uy.^5.*dxux.*dtT + 8.*T.*ux.^2.*uy.^5.*dxux.*dtT + 4.*T.*uy.^7.*dxux.*dtT - 8.*T.*ux.*dxuy.*dtT - 16.*T.*ux.^3.*dxuy.*dtT - 10.*T.*ux.^5.*dxuy.*dtT - 2.*T.*ux.^7.*dxuy.*dtT - 24.*T.*ux.*uy.^2.*dxuy.*dtT - 32.*T.*ux.^3.*uy.^2.*dxuy.*dtT - 10.*T.*ux.^5.*uy.^2.*dxuy.*dtT - 22.*T.*ux.*uy.^4.*dxuy.*dtT - 14.*T.*ux.^3.*uy.^4.*dxuy.*dtT - 6.*T.*ux.*uy.^6.*dxuy.*dtT + 16.*pi.*T.^3.*ux.*uy.*dtux + 32.*pi.*T.^3.*ux.^3.*uy.*dtux + 20.*pi.*T.^3.*ux.^5.*uy.*dtux + 4.*pi.*T.^3.*ux.^7.*uy.*dtux + 32.*pi.*T.^3.*ux.*uy.^3.*dtux + 40.*pi.*T.^3.*ux.^3.*uy.^3.*dtux + 12.*pi.*T.^3.*ux.^5.*uy.^3.*dtux + 20.*pi.*T.^3.*ux.*uy.^5.*dtux + 12.*pi.*T.^3.*ux.^3.*uy.^5.*dtux + 4.*pi.*T.^3.*ux.*uy.^7.*dtux + 8.*T.*ux.*dyT.*dtux + 12.*T.*ux.^3.*dyT.*dtux + 6.*T.*ux.^5.*dyT.*dtux + T.*ux.^7.*dyT.*dtux + 16.*T.*ux.*uy.^2.*dyT.*dtux + 14.*T.*ux.^3.*uy.^2.*dyT.*dtux + 2.*T.*ux.^5.*uy.^2.*dyT.*dtux + 8.*T.*ux.*uy.^4.*dyT.*dtux + T.*ux.^3.*uy.^4.*dyT.*dtux - 4.*T.^2.*dyux.*dtux - 4.*T.^2.*ux.^2.*dyux.*dtux - T.^2.*ux.^4.*dyux.*dtux - 12.*T.^2.*uy.^2.*dyux.*dtux - 4.*T.^2.*ux.^2.*uy.^2.*dyux.*dtux + 3.*T.^2.*ux.^4.*uy.^2.*dyux.*dtux - 11.*T.^2.*uy.^4.*dyux.*dtux - 3.*T.^2.*uy.^6.*dyux.*dtux + 8.*T.^2.*ux.*uy.*dyuy.*dtux + 6.*T.^2.*ux.^3.*uy.*dyuy.*dtux + 14.*T.^2.*ux.*uy.^3.*dyuy.*dtux + 6.*T.^2.*ux.^3.*uy.^3.*dyuy.*dtux + 6.*T.^2.*ux.*uy.^5.*dyuy.*dtux + 6.*T.*ux.^2.*uy.*dxT.*dtux + 5.*T.*ux.^4.*uy.*dxT.*dtux + 2.*T.*uy.^3.*dxT.*dtux + 8.*T.*ux.^2.*uy.^3.*dxT.*dtux + T.*ux.^4.*uy.^3.*dxT.*dtux + 3.*T.*uy.^5.*dxT.*dtux + 2.*T.*ux.^2.*uy.^5.*dxT.*dtux + T.*uy.^7.*dxT.*dtux - 8.*T.^2.*ux.*uy.*dxux.*dtux - 6.*T.^2.*ux.^3.*uy.*dxux.*dtux - 14.*T.^2.*ux.*uy.^3.*dxux.*dtux - 6.*T.^2.*ux.^3.*uy.^3.*dxux.*dtux - 6.*T.^2.*ux.*uy.^5.*dxux.*dtux - 4.*T.^2.*dxuy.*dtux - 4.*T.^2.*ux.^2.*dxuy.*dtux - T.^2.*ux.^4.*dxuy.*dtux - 12.*T.^2.*uy.^2.*dxuy.*dtux - 4.*T.^2.*ux.^2.*uy.^2.*dxuy.*dtux + 3.*T.^2.*ux.^4.*uy.^2.*dxuy.*dtux - 11.*T.^2.*uy.^4.*dxuy.*dtux - 3.*T.^2.*uy.^6.*dxuy.*dtux + 16.*pi.*T.^3.*dtuy + 48.*pi.*T.^3.*ux.^2.*dtuy + 52.*pi.*T.^3.*ux.^4.*dtuy + 24.*pi.*T.^3.*ux.^6.*dtuy + 4.*pi.*T.^3.*ux.^8.*dtuy + 64.*pi.*T.^3.*uy.^2.*dtuy + 136.*pi.*T.^3.*ux.^2.*uy.^2.*dtuy + 92.*pi.*T.^3.*ux.^4.*uy.^2.*dtuy + 20.*pi.*T.^3.*ux.^6.*uy.^2.*dtuy + 84.*pi.*T.^3.*uy.^4.*dtuy + 112.*pi.*T.^3.*ux.^2.*uy.^4.*dtuy + 36.*pi.*T.^3.*ux.^4.*uy.^4.*dtuy + 44.*pi.*T.^3.*uy.^6.*dtuy + 28.*pi.*T.^3.*ux.^2.*uy.^6.*dtuy + 8.*pi.*T.^3.*uy.^8.*dtuy + 8.*T.*uy.*dyT.*dtuy + 16.*T.*ux.^2.*uy.*dyT.*dtuy + 12.*T.*ux.^4.*uy.*dyT.*dtuy + 3.*T.*ux.^6.*uy.*dyT.*dtuy + 20.*T.*uy.^3.*dyT.*dtuy + 26.*T.*ux.^2.*uy.^3.*dyT.*dtuy + 8.*T.*ux.^4.*uy.^3.*dyT.*dtuy + 14.*T.*uy.^5.*dyT.*dtuy + 7.*T.*ux.^2.*uy.^5.*dyT.*dtuy + 2.*T.*uy.^7.*dyT.*dtuy - 4.*T.^2.*ux.*uy.*dyux.*dtuy - 8.*T.^2.*ux.^3.*uy.*dyux.*dtuy - 3.*T.^2.*ux.^5.*uy.*dyux.*dtuy + 3.*T.^2.*ux.*uy.^5.*dyux.*dtuy - 4.*T.^2.*dyuy.*dtuy - 6.*T.^2.*ux.^2.*dyuy.*dtuy - 2.*T.^2.*ux.^4.*dyuy.*dtuy - 10.*T.^2.*uy.^2.*dyuy.*dtuy - 16.*T.^2.*ux.^2.*uy.^2.*dyuy.*dtuy - 6.*T.^2.*ux.^4.*uy.^2.*dyuy.*dtuy - 6.*T.^2.*uy.^4.*dyuy.*dtuy - 6.*T.^2.*ux.^2.*uy.^4.*dyuy.*dtuy + 2.*T.*ux.^3.*dxT.*dtuy + 3.*T.*ux.^5.*dxT.*dtuy + T.*ux.^7.*dxT.*dtuy + 6.*T.*ux.*uy.^2.*dxT.*dtuy + 8.*T.*ux.^3.*uy.^2.*dxT.*dtuy + 2.*T.*ux.^5.*uy.^2.*dxT.*dtuy + 5.*T.*ux.*uy.^4.*dxT.*dtuy + T.*ux.^3.*uy.^4.*dxT.*dtuy + 4.*T.^2.*dxux.*dtuy + 6.*T.^2.*ux.^2.*dxux.*dtuy + 2.*T.^2.*ux.^4.*dxux.*dtuy + 10.*T.^2.*uy.^2.*dxux.*dtuy + 16.*T.^2.*ux.^2.*uy.^2.*dxux.*dtuy + 6.*T.^2.*ux.^4.*uy.^2.*dxux.*dtuy + 6.*T.^2.*uy.^4.*dxux.*dtuy + 6.*T.^2.*ux.^2.*uy.^4.*dxux.*dtuy - 4.*T.^2.*ux.*uy.*dxuy.*dtuy - 8.*T.^2.*ux.^3.*uy.*dxuy.*dtuy - 3.*T.^2.*ux.^5.*uy.*dxuy.*dtuy + 3.*T.^2.*ux.*uy.^5.*dxuy.*dtuy + 4.*T.*ux.^2.*dydtT + 8.*T.*ux.^4.*dydtT + 5.*T.*ux.^6.*dydtT + T.*ux.^8.*dydtT + 4.*T.*uy.^2.*dydtT + 18.*T.*ux.^2.*uy.^2.*dydtT + 18.*T.*ux.^4.*uy.^2.*dydtT + 5.*T.*ux.^6.*uy.^2.*dydtT + 10.*T.*uy.^4.*dydtT + 21.*T.*ux.^2.*uy.^4.*dydtT + 9.*T.*ux.^4.*uy.^4.*dydtT + 8.*T.*uy.^6.*dydtT + 7.*T.*ux.^2.*uy.^6.*dydtT + 2.*T.*uy.^8.*dydtT - 4.*T.^2.*ux.*dydtux - 8.*T.^2.*ux.^3.*dydtux - 5.*T.^2.*ux.^5.*dydtux - T.^2.*ux.^7.*dydtux - 12.*T.^2.*ux.*uy.^2.*dydtux - 16.*T.^2.*ux.^3.*uy.^2.*dydtux - 5.*T.^2.*ux.^5.*uy.^2.*dydtux - 11.*T.^2.*ux.*uy.^4.*dydtux - 7.*T.^2.*ux.^3.*uy.^4.*dydtux - 3.*T.^2.*ux.*uy.^6.*dydtux - 4.*T.^2.*uy.*dydtuy - 6.*T.^2.*ux.^2.*uy.*dydtuy - 2.*T.^2.*ux.^4.*uy.*dydtuy - 10.*T.^2.*uy.^3.*dydtuy - 10.*T.^2.*ux.^2.*uy.^3.*dydtuy - 2.*T.^2.*ux.^4.*uy.^3.*dydtuy - 8.*T.^2.*uy.^5.*dydtuy - 4.*T.^2.*ux.^2.*uy.^5.*dydtuy - 2.*T.^2.*uy.^7.*dydtuy + 2.*T.*ux.^3.*uy.*dxdtT + 3.*T.*ux.^5.*uy.*dxdtT + T.*ux.^7.*uy.*dxdtT + 2.*T.*ux.*uy.^3.*dxdtT + 6.*T.*ux.^3.*uy.^3.*dxdtT + 3.*T.*ux.^5.*uy.^3.*dxdtT + 3.*T.*ux.*uy.^5.*dxdtT + 3.*T.*ux.^3.*uy.^5.*dxdtT + T.*ux.*uy.^7.*dxdtT + 4.*T.^2.*uy.*dxdtux + 6.*T.^2.*ux.^2.*uy.*dxdtux + 2.*T.^2.*ux.^4.*uy.*dxdtux + 10.*T.^2.*uy.^3.*dxdtux + 10.*T.^2.*ux.^2.*uy.^3.*dxdtux + 2.*T.^2.*ux.^4.*uy.^3.*dxdtux + 8.*T.^2.*uy.^5.*dxdtux + 4.*T.^2.*ux.^2.*uy.^5.*dxdtux + 2.*T.^2.*uy.^7.*dxdtux - 4.*T.^2.*ux.*dxdtuy - 8.*T.^2.*ux.^3.*dxdtuy - 5.*T.^2.*ux.^5.*dxdtuy - T.^2.*ux.^7.*dxdtuy - 12.*T.^2.*ux.*uy.^2.*dxdtuy - 16.*T.^2.*ux.^3.*uy.^2.*dxdtuy - 5.*T.^2.*ux.^5.*uy.^2.*dxdtuy - 11.*T.^2.*ux.*uy.^4.*dxdtuy - 7.*T.^2.*ux.^3.*uy.^4.*dxdtuy - 3.*T.^2.*ux.*uy.^6.*dxdtuy));

    res = [L1 L2 L3] + R;

    % compute the gradients. hyhyhy



    dTL1_t1 = (1./(9.*(2 + ux.^2 + uy.^2))).*(16.* pi.^2.*(8.*pi.*T.^2 + 16.*pi.*T.^2.*ux.^2 + 6.*pi.*T.^2.*ux.^4 + 16.*pi.*T.^2.*uy.^2 + 12.*pi.*T.^2.*ux.^2.*uy.^2 + 6.*pi.*T.^2.*uy.^4 + ux.^2.*uy.*dyT + uy.^3.*dyT - 4.*T.*ux.*uy.* dyux + 2.*T.*ux.^2.* dyuy - 2.*T.*uy.^2.*dyuy + ux.^3.*dxT + ux.*uy.^2.*dxT - 2.*T.*ux.^2.*dxux + 2.*T.*uy.^2.* dxux - 4.*T.*ux.*uy.* dxuy));
    dTL1_t2 = (16.*pi.^2.*T.*ux.*(ux.^2 + uy.^2))./(9.*(2 + ux.^2 + uy.^2));
    dTL1_t3 = (16.*pi.^2.*T.*uy.*(ux.^2 + uy.^2))./(9.*(2 + ux.^2 + uy.^2));

    dTL2_t1 = (8.*pi.^2.*(24.*pi.*T.^2.*ux + 36.*pi.*T.^2.*ux.^3 + 12.*pi.*T.^2.*ux.^5 + 36.*pi.*T.^2.*ux.*uy.^2 + 24.*pi.*T.^2.*ux.^3.*uy.^2 + 12.*pi.*T.^2.*ux.*uy.^4 + ux.^3.*uy.*dyT + ux.*uy.^3.*dyT - 4.*T.*uy.*dyux - 6.*T.*ux.^2.*uy .* dyux - 2.*T.*uy.^3.*dyux + 4.*T.*ux.*dyuy + 4.*T.*ux.^3.*dyuy + 2.*ux.^2.*dxT + 2.*ux.^4.*dxT + 2.*uy.^2.*dxT + 3.*ux.^2.*uy.^2.*dxT + uy.^4.*dxT - 4.*T.*ux.*dxux - 4.*T.*ux.^3.*dxux - 4.*T.*uy.*dxuy - 6.*T.*ux.^2.*uy.* dxuy - 2.*T.*uy.^3.*dxuy))./(9.*sqrt(1 + ux.^2 + uy.^2).*(2 + ux.^2 + uy.^2));
    dTL2_t2 = (8.*pi.^2.*T.*(ux.^2 + uy.^2).*(2 + 2.*ux.^2 + uy.^2))./ (9.*sqrt(1 + ux.^2 + uy.^2).*(2 + ux.^2 + uy.^2));
    dTL2_t3 = (8.*pi.^2.*T.*ux.* uy.*(ux.^2 + uy.^2))./ (9.*sqrt(1 + ux.^2 + uy.^2).*(2 + ux.^2 + uy.^2));
 
    dTL3_t1 = (8.*pi.^2.*(24.*pi.*T.^2.*uy + 36.*pi.*T.^2.*ux.^2.*uy + 12.*pi.*T.^2.*ux.^4.*uy + 36.*pi.*T.^2.*uy.^3 + 24.*pi.*T.^2.*ux.^2.*uy.^3 + 12.*pi.*T.^2.*uy.^5 + 2.*ux.^2.*dyT + ux.^4.*dyT + 2.*uy.^2.*dyT + 3.*ux.^2.*uy.^2.*dyT + 2.*uy.^4.*dyT - 4.*T.*ux.*dyux - 2.*T.*ux.^3.*dyux - 6.*T.*ux.*uy.^2.*dyux - 4.*T.*uy.*dyuy - 4.*T.*uy.^3.*dyuy + ux.^3.*uy.*dxT + ux.*uy.^3.*dxT + 4.*T.*uy.*dxux + 4.*T.*uy.^3.*dxux - 4.*T.*ux.*dxuy - 2.*T.*ux.^3.*dxuy - 6.*T.*ux.*uy.^2.*dxuy))./(9.*sqrt(1 + ux.^2 + uy.^2).*(2 + ux.^2 + uy.^2));
    dTL3_t2 = (8.*pi.^2.*T.*ux.*uy.*(ux.^2 + uy.^2))./ (9.*sqrt(1 + ux.^2 + uy.^2).*(2 + ux.^2 + uy.^2));
    dTL3_t3 = (8.*pi.^2.*T.*(ux.^2 + uy.^2).*(2 + ux.^2 + 2.*uy.^2))./ (9.*sqrt(1 + ux.^2 + uy.^2).*(2 + ux.^2 + uy.^2));

    % ********************

    duxL1_t1 = (1./(9.*(2 + ux.^2 + uy.^2).^2)).*(16.*pi.^2.* T.*(16.*pi.*T.^2.*ux + 16.*pi.*T.^2.*ux.^3 + 4.*pi.*T.^2.*ux.^5 + 16.*pi.*T.^2.*ux.*uy.^2 + 8.*pi.*T.^2.*ux.^3.*uy.^2 + 4.*pi.*T.^2.*ux.*uy.^4 + 4.*ux.*uy.*dyT - 4.*T.*uy.*dyux + 2.*T.*ux.^2.*uy.* dyux - 2.*T.*uy.^3.*dyux + 4.*T.*ux.*dyuy + 4.*T.*ux.*uy.^2.* dyuy + 6.*ux.^2.*dxT + ux.^4.*dxT + 2.*uy.^2.*dxT + 2.*ux.^2.*uy.^2.* dxT + uy.^4.*dxT - 4.*T.*ux.*dxux - 4.*T.*ux.*uy.^2.* dxux - 4.*T.*uy.*dxuy + 2.*T.*ux.^2.*uy.* dxuy - 2.*T.*uy.^3.*dxuy));
    duxL1_t2 = (1./(9.*(2 + ux.^2 + uy.^2).^2)).*(16.*pi.^2.* T.*(- 2.*T.*ux.^2 - T.*ux.^4 + 2.*T.*uy.^2 + T.*uy.^4));
    duxL1_t3 = (1./(9.*(2 + ux.^2 + uy.^2).^2)).*(16.*pi.^2.* T.*(-4.*T.*ux.*uy - 2.*T.*ux.^3.*uy - 2.*T.*ux.*uy.^3));

    duxL2_t1 = (8.*pi.^2.*T.*(16.*pi.*T.^2 + 64.*pi.*T.^2.*ux.^2 + 84.*pi.*T.^2.*ux.^4 + 44.*pi.*T.^2.*ux.^6 + 8.*pi.*T.^2.*ux.^8 + 48.*pi.*T.^2.*uy.^2 + 136.*pi.*T.^2.*ux.^2.*uy.^2 + 112.*pi.*T.^2.*ux.^4.*uy.^2 + 28.*pi.*T.^2.*ux.^6.*uy.^2 + 52.*pi.*T.^2.*uy.^4 + 92.*pi.*T.^2.*ux.^2.*uy.^4 + 36.*pi.*T.^2.*ux.^4.*uy.^4 + 24.*pi.*T.^2.*uy.^6 + 20.*pi.*T.^2.*ux.^2.*uy.^6 + 4.*pi.*T.^2.*uy.^8 + 6.*ux.^2.*uy.* dyT + 5.*ux.^4.*uy.*dyT + 2.*uy.^3.*dyT + 8.*ux.^2.*uy.^3.* dyT + ux.^4.*uy.^3.*dyT + 3.*uy.^5.*dyT + 2.*ux.^2.*uy.^5.* dyT + uy.^7.*dyT - 4.*T.*ux.*uy.* dyux + 3.*T.*ux.^5.*uy.* dyux - 8.*T.*ux.*uy.^3.* dyux - 3.*T.*ux.*uy.^5.* dyux + 4.*T.*dyuy + 10.*T.*ux.^2.*dyuy + 6.*T.*ux.^4.*dyuy + 6.*T.*uy.^2.* dyuy + 16.*T.*ux.^2.*uy.^2.* dyuy + 6.*T.*ux.^4.*uy.^2.* dyuy + 2.*T.*uy.^4.*dyuy + 6.*T.*ux.^2.*uy.^4.* dyuy + 8.*ux.*dxT + 20.*ux.^3.*dxT + 14.*ux.^5.*dxT + 2.*ux.^7.*dxT + 16.*ux.*uy.^2.*dxT + 26.*ux.^3.*uy.^2.* dxT + 7.*ux.^5.*uy.^2.*dxT + 12.*ux.*uy.^4.* dxT + 8.*ux.^3.*uy.^4.*dxT + 3.*ux.*uy.^6.* dxT - 4.*T.*dxux - 10.*T.*ux.^2.*dxux - 6.*T.*ux.^4.* dxux - 6.*T.*uy.^2.*dxux - 16.*T.*ux.^2.*uy.^2.* dxux - 6.*T.*ux.^4.*uy.^2.* dxux - 2.*T.*uy.^4.* dxux - 6.*T.*ux.^2.*uy.^4.* dxux - 4.*T.*ux.*uy.* dxuy + 3.*T.*ux.^5.*uy.* dxuy - 8.*T.*ux.*uy.^3.* dxuy - 3.*T.*ux.*uy.^5.* dxuy))./(9.*(1 + ux.^2 + uy.^2).^(3./2).*(2 + ux.^2 + uy.^2).^2);
    duxL2_t2 = -((16.*pi.^2.*T.^2.*ux.*(1 + ux.^2))./(9.*sqrt(1 + ux.^2 + uy.^2).*(2 + ux.^2 + uy.^2)));
    duxL2_t3 = -((8.*pi.^2.*T.^2.* uy.*(2 + 3.*ux.^2 + uy.^2))./(9.* sqrt(1 + ux.^2 + uy.^2).* (2 + ux.^2 + uy.^2)));

    duxL3_t1 = (8.*pi.^2.*T.*(16.*pi.*T.^2.*ux.*uy + 32.*pi.*T.^2.*ux.^3.*uy + 20.*pi.*T.^2.*ux.^5.*uy + 4.*pi.*T.^2.*ux.^7.*uy + 32.*pi.*T.^2.*ux.*uy.^3 + 40.*pi.*T.^2.*ux.^3.*uy.^3 + 12.*pi.*T.^2.*ux.^5.*uy.^3 + 20.*pi.*T.^2.*ux.*uy.^5 + 12.*pi.*T.^2.*ux.^3.*uy.^5 + 4.*pi.*T.^2.*ux.*uy.^7 + 8.*ux.*dyT + 12.*ux.^3.*dyT + 6.*ux.^5.*dyT + ux.^7.*dyT + 16.*ux.*uy.^2.*dyT + 14.*ux.^3.*uy.^2.* dyT + 2.*ux.^5.*uy.^2.*dyT + 8.*ux.*uy.^4.* dyT + ux.^3.*uy.^4.*dyT - 4.*T.*dyux - 4.*T.*ux.^2.* dyux - T.*ux.^4.*dyux - 12.*T.*uy.^2.* dyux - 4.*T.*ux.^2.*uy.^2.* dyux + 3.*T.*ux.^4.*uy.^2.* dyux - 11.*T.*uy.^4.*dyux - 3.*T.*uy.^6.* dyux + 8.*T.*ux.*uy.* dyuy + 6.*T.*ux.^3.*uy.* dyuy + 14.*T.*ux.*uy.^3.* dyuy + 6.*T.*ux.^3.*uy.^3.* dyuy + 6.*T.*ux.*uy.^5.* dyuy + 6.*ux.^2.*uy.* dxT + 5.*ux.^4.*uy.*dxT + 2.*uy.^3.*dxT + 8.*ux.^2.*uy.^3.* dxT + ux.^4.*uy.^3.*dxT + 3.*uy.^5.*dxT + 2.*ux.^2.*uy.^5.* dxT + uy.^7.*dxT - 8.*T.*ux.*uy.* dxux - 6.*T.*ux.^3.*uy.* dxux - 14.*T.*ux.*uy.^3.* dxux - 6.*T.*ux.^3.*uy.^3.* dxux - 6.*T.*ux.*uy.^5.* dxux - 4.*T.*dxuy - 4.*T.*ux.^2.* dxuy - T.*ux.^4.*dxuy - 12.*T.*uy.^2.* dxuy - 4.*T.*ux.^2.*uy.^2.* dxuy + 3.*T.*ux.^4.*uy.^2.* dxuy - 11.*T.*uy.^4.*dxuy - 3.*T.*uy.^6.*dxuy))./(9.*(1 + ux.^2 + uy.^2).^(3./2).*(2 + ux.^2 + uy.^2).^2);
    duxL3_t2 = (16.*pi.^2.*T.^2.*uy.*(1 + uy.^2))./(9.*sqrt(1 + ux.^2 + uy.^2).*(2 + ux.^2 + uy.^2));
    duxL3_t3 = -((8.*pi.^2.*T.^2.* ux.*(2 + ux.^2 + 3.*uy.^2))./(9.* sqrt(1 + ux.^2 + uy.^2).* (2 + ux.^2 + uy.^2)));

    % ********************
    
    duyL1_t1 = (1./(9.*(2 + ux.^2 + uy.^2).^2)).*(16.*pi.^2.* T.*(16.*pi.*T.^2.*uy + 16.*pi.*T.^2.*ux.^2.*uy + 4.*pi.*T.^2.*ux.^4.*uy + 16.*pi.*T.^2.*uy.^3 + 8.*pi.*T.^2.*ux.^2.*uy.^3 + 4.*pi.*T.^2.*uy.^5 + 2.*ux.^2.*dyT + ux.^4.*dyT + 6.*uy.^2.*dyT + 2.*ux.^2.*uy.^2.* dyT + uy.^4.*dyT - 4.*T.*ux.*dyux - 2.*T.*ux.^3.* dyux + 2.*T.*ux.*uy.^2.* dyux - 4.*T.*uy.*dyuy - 4.*T.*ux.^2.*uy.* dyuy + 4.*ux.*uy.*dxT + 4.*T.*uy.*dxux + 4.*T.*ux.^2.*uy.* dxux - 4.*T.*ux.*dxuy - 2.*T.*ux.^3.* dxuy + 2.*T.*ux.*uy.^2.* dxuy));
    duyL1_t2 = -((32.*pi.^2.*T.^2.*ux.*uy)./(9.*(2 + ux.^2 + uy.^2)));
    duyL1_t3 = (16.*pi.^2.*T.^2.*(ux - uy).*(ux + uy))./(9.*(2 + ux.^2 + uy.^2));

    duyL2_t1 = (8.*pi.^2.*T.*(16.*pi.*T.^2.*ux.*uy + 32.*pi.*T.^2.*ux.^3.*uy + 20.*pi.*T.^2.*ux.^5.*uy + 4.*pi.*T.^2.*ux.^7.*uy + 32.*pi.*T.^2.*ux.*uy.^3 + 40.*pi.*T.^2.*ux.^3.*uy.^3 + 12.*pi.*T.^2.*ux.^5.*uy.^3 + 20.*pi.*T.^2.*ux.*uy.^5 + 12.*pi.*T.^2.*ux.^3.*uy.^5 + 4.*pi.*T.^2.*ux.*uy.^7 + 2.*ux.^3.*dyT + 3.*ux.^5.*dyT + ux.^7.*dyT + 6.*ux.*uy.^2.*dyT + 8.*ux.^3.*uy.^2.*dyT + 2.*ux.^5.*uy.^2.* dyT + 5.*ux.*uy.^4.*dyT + ux.^3.*uy.^4.* dyT - 4.*T.*dyux - 12.*T.*ux.^2.*dyux - 11.*T.*ux.^4.* dyux - 3.*T.*ux.^6.*dyux - 4.*T.*uy.^2.* dyux - 4.*T.*ux.^2.*uy.^2.* dyux - T.*uy.^4.*dyux + 3.*T.*ux.^2.*uy.^4.* dyux - 8.*T.*ux.*uy.* dyuy - 14.*T.*ux.^3.*uy.* dyuy - 6.*T.*ux.^5.*uy.* dyuy - 6.*T.*ux.*uy.^3.* dyuy - 6.*T.*ux.^3.*uy.^3.* dyuy + 8.*uy.*dxT + 16.*ux.^2.*uy.* dxT + 8.*ux.^4.*uy.*dxT + 12.*uy.^3.*dxT + 14.*ux.^2.*uy.^3.* dxT + ux.^4.*uy.^3.*dxT + 6.*uy.^5.*dxT + 2.*ux.^2.*uy.^5.* dxT + uy.^7.*dxT + 8.*T.*ux.*uy.* dxux + 14.*T.*ux.^3.*uy.* dxux + 6.*T.*ux.^5.*uy.* dxux + 6.*T.*ux.*uy.^3.* dxux + 6.*T.*ux.^3.*uy.^3.* dxux - 4.*T.*dxuy - 12.*T.*ux.^2.* dxuy - 11.*T.*ux.^4.*dxuy - 3.*T.*ux.^6.* dxuy - 4.*T.*uy.^2.*dxuy - 4.*T.*ux.^2.*uy.^2.* dxuy - T.*uy.^4.*dxuy + 3.*T.*ux.^2.*uy.^4.* dxuy))./ (9.*(1 + ux.^2 + uy.^2).^(3./2).*(2 + ux.^2 + uy.^2).^2);
    duyL2_t2 = -((8.*pi.^2.*T.^2.* uy.*(2 + 3.*ux.^2 + uy.^2))./(9.* sqrt(1 + ux.^2 + uy.^2).* (2 + ux.^2 + uy.^2)));
    duyL2_t3 = (16.*pi.^2.*T.^2.*ux.*(1 + ux.^2))./(9.* sqrt(1 + ux.^2 + uy.^2).*(2 + ux.^2 + uy.^2));


    duyL3_t1 = (8.*pi.^2.*T.*(16.*pi.*T.^2 + 48.*pi.*T.^2.*ux.^2 + 52.*pi.*T.^2.*ux.^4 + 24.*pi.*T.^2.*ux.^6 + 4.*pi.*T.^2.*ux.^8 + 64.*pi.*T.^2.*uy.^2 + 136.*pi.*T.^2.*ux.^2.*uy.^2 + 92.*pi.*T.^2.*ux.^4.*uy.^2 + 20.*pi.*T.^2.*ux.^6.*uy.^2 + 84.*pi.*T.^2.*uy.^4 + 112.*pi.*T.^2.*ux.^2.*uy.^4 + 36.*pi.*T.^2.*ux.^4.*uy.^4 + 44.*pi.*T.^2.*uy.^6 + 28.*pi.*T.^2.*ux.^2.*uy.^6 + 8.*pi.*T.^2.*uy.^8 + 8.*uy.*dyT + 16.*ux.^2.*uy.*dyT + 12.*ux.^4.*uy.*dyT + 3.*ux.^6.*uy.* dyT + 20.*uy.^3.*dyT + 26.*ux.^2.*uy.^3.* dyT + 8.*ux.^4.*uy.^3.* dyT + 14.*uy.^5.*dyT + 7.*ux.^2.*uy.^5.*dyT + 2.*uy.^7.*dyT - 4.*T.*ux.*uy.* dyux - 8.*T.*ux.^3.*uy.* dyux - 3.*T.*ux.^5.*uy.* dyux + 3.*T.*ux.*uy.^5.* dyux - 4.*T.*dyuy - 6.*T.*ux.^2.* dyuy - 2.*T.*ux.^4.*dyuy - 10.*T.*uy.^2.* dyuy - 16.*T.*ux.^2.*uy.^2.* dyuy - 6.*T.*ux.^4.*uy.^2.* dyuy - 6.*T.*uy.^4.*dyuy - 6.*T.*ux.^2.*uy.^4.* dyuy + 2.*ux.^3.*dxT + 3.*ux.^5.*dxT + ux.^7.*dxT + 6.*ux.*uy.^2.*dxT + 8.*ux.^3.*uy.^2.*dxT + 2.*ux.^5.*uy.^2.* dxT + 5.*ux.*uy.^4.*dxT + ux.^3.*uy.^4.* dxT + 4.*T.*dxux + 6.*T.*ux.^2.*dxux + 2.*T.*ux.^4.* dxux + 10.*T.*uy.^2.*dxux + 16.*T.*ux.^2.*uy.^2.* dxux + 6.*T.*ux.^4.*uy.^2.* dxux + 6.*T.*uy.^4.* dxux + 6.*T.*ux.^2.*uy.^4.* dxux - 4.*T.*ux.*uy.* dxuy - 8.*T.*ux.^3.*uy.* dxuy - 3.*T.*ux.^5.*uy.* dxuy + 3.*T.*ux.*uy.^5.* dxuy))./(9.*(1 + ux.^2 + uy.^2).^(3./2).*(2 + ux.^2 + uy.^2).^2);
    duyL3_t2 = -((8.*pi.^2.*T.^2.* ux.*(2 + ux.^2 + 3.*uy.^2))./(9.* sqrt(1 + ux.^2 + uy.^2).* (2 + ux.^2 + uy.^2)));
    duyL3_t3 = -((16.*pi.^2.*T.^2.*uy.*(1 + uy.^2))./(9.* sqrt(1 + ux.^2 + uy.^2).*(2 + ux.^2 + uy.^2)));

    J = zeros(3*GRID_SIZE*GRID_SIZE, 3*GRID_SIZE*GRID_SIZE);
    for i=1:GRID_SIZE*GRID_SIZE
      x = rem(i, GRID_SIZE);
      y = mod(i, GRID_SIZE);

      du = sparse(zeros(GRID_SIZE, GRID_SIZE));
      du(i) = 1;
      dxdu = dx(du);
      dydu = dy(du);

      Jxy = dTL1_t1 .* du + dTL1_t2 .* dxdu + dTL1_t3 .* dydu;
      J(                      1:  GRID_SIZE*GRID_SIZE,i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dTL2_t1 .* du + dTL2_t2 .* dxdu + dTL2_t3 .* dydu;
      J(  GRID_SIZE*GRID_SIZE+1:2*GRID_SIZE*GRID_SIZE,i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = dTL3_t1 .* du + dTL3_t2 .* dxdu + dTL3_t3 .* dydu;
      J(2*GRID_SIZE*GRID_SIZE+1:3*GRID_SIZE*GRID_SIZE,i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);

      Jxy = duxL1_t1 .* du + duxL1_t2 .* dxdu + duxL1_t3 .* dydu;
      J(                      1:  GRID_SIZE*GRID_SIZE,  GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = duxL2_t1 .* du + duxL2_t2 .* dxdu + duxL2_t3 .* dydu;
      J(  GRID_SIZE*GRID_SIZE+1:2*GRID_SIZE*GRID_SIZE,  GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = duxL3_t1 .* du + duxL3_t2 .* dxdu + duxL3_t3 .* dydu;
      J(2*GRID_SIZE*GRID_SIZE+1:3*GRID_SIZE*GRID_SIZE,  GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);

      Jxy = duyL1_t1 .* du + duyL1_t2 .* dxdu + duyL1_t3 .* dydu;
      J(                      1:  GRID_SIZE*GRID_SIZE,2*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = duyL2_t1 .* du + duyL2_t2 .* dxdu + duyL2_t3 .* dydu;
      J(  GRID_SIZE*GRID_SIZE+1:2*GRID_SIZE*GRID_SIZE,2*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
      Jxy = duyL3_t1 .* du + duyL3_t2 .* dxdu + duyL3_t3 .* dydu;
      J(2*GRID_SIZE*GRID_SIZE+1:3*GRID_SIZE*GRID_SIZE,2*GRID_SIZE*GRID_SIZE+i) = reshape(Jxy, GRID_SIZE*GRID_SIZE, 1);
    end
    % sqrt(sum(sum(res.*res)))
  end

  function r = phydro()
    frame_active = (mod(t-1,int16(1/DT/SECONDS_PER_LIGHTYEAR/30)) == 0);
    if VIDEO & ~frame_active 
      return
    end
    if VIDEO
      disp(t);
    end
    plotsx = 4;
    plotsy = 3;

    subplot(plotsy, plotsx, 1);
    image((T-TEMP_OFFSET) * 256 * 100 + 128); 
    
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vt,frame);
    % end

    subplot(plotsy, plotsx, 2);
    image(ux*128 + 128);
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vvx,frame);
    % end

    subplot(plotsy, plotsx, 3);
    image(uy*128 + 128);
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vvy,frame);
    % end

    subplot(plotsy, plotsx, 4);
    image(sqrt(ux.*ux + uy.*uy)*128 + 128);
    % if (VIDEO) & (frame_active)
    %   frame = getframe;
    %   writeVideo(vv,frame);
    % end

    subplot(plotsy, plotsx, 5);
    image(dtT*4+128);

    subplot(plotsy, plotsx, 6);
    image(dtux*4+128);
    subplot(plotsy, plotsx, 7);
    image(dtuy*4+128);

    %truesize(fig);
    if (VIDEO) & (frame_active)
      frame = getframe(gcf);
      writeVideo(vt,frame);
    end

    
    % subplot(plotsy, plotsx, 9);
    % image(dtdtvx*8 + 128);
    % subplot(plotsy, plotsx, 10);
    % image(dtdtvy*8 + 128);
    % subplot(plotsy, plotsx, 11);
    % image(Pi*2 + 128);
    % subplot(plotsy, plotsx, 12);
    % image(B*2 + 128);
    return;

  end

  % [-1 0 1
  function r = dx2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for y = 1:GRID_SIZE
      r(1,y) = - arg(GRID_SIZE,y)/2 + arg(2, y)/2;
    end
    for y = 1:GRID_SIZE
      r(GRID_SIZE,y) =  -arg(GRID_SIZE-1,y)/2 + arg(1, y)/2;
    end
    for x = 2:GRID_SIZE-1
      for y = 1:GRID_SIZE
        r(x, y) = (arg(x+1,y) - arg(x-1,y))/2;   
      end
    end 
    r = r / DX;
  end
  function r = dy2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for x = 1:GRID_SIZE
      r(x,1) = - arg(x,GRID_SIZE)/2 + arg(x,2)/2;
    end
    for x = 1:GRID_SIZE
      r(x,GRID_SIZE) =  -arg(x,GRID_SIZE-1)/2 + arg(x,1)/2;
    end
    for x = 1:GRID_SIZE
      for y = 2:GRID_SIZE-1
        r(x, y) = (arg(x,y+1) - arg(x,y-1))/2;   
      end
    end 
    r = r / DX;
  end
  function r = dx(arg)
    r = d1 * arg;
  end
  function r = dy(arg)
    r = arg * (d1.');
  end
  function r = dxdx(arg)
    r = d2 * arg;
  end
  function r = dydy(arg)
    r = arg * (d2.');
  end
  function r = dxdx2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for y = 1:GRID_SIZE
      r(1,y) = arg(GRID_SIZE,y) + arg(2, y) - 2*arg(1,y);
    end
    for y = 1:GRID_SIZE
      r(GRID_SIZE,y) =  arg(GRID_SIZE-1,y) + arg(1, y) - 2*arg(GRID_SIZE,y);
    end
    for x = 2:GRID_SIZE-1
      for y = 1:GRID_SIZE
        r(x, y) = arg(x+1,y) + arg(x-1,y) - 2*arg(x,y);
      end
    end 
    r = r / DX / DX;
  end
  function r = dydy2(arg)
    r = zeros(GRID_SIZE, GRID_SIZE);
    for x = 1:GRID_SIZE
      r(x,1) = arg(x,GRID_SIZE) + arg(x,2) - 2*arg(x,1);
    end
    for x = 1:GRID_SIZE
      r(x,GRID_SIZE) =  arg(x,GRID_SIZE-1) + arg(x,1) - 2*arg(x,GRID_SIZE);
    end
    for x = 1:GRID_SIZE
      for y = 2:GRID_SIZE-1
        r(x, y) = arg(x,y+1) + arg(x,y-1) - 2*arg(x,y);   
      end
    end 
    r = r / DX / DX;
  end

end


